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ABSTRACT 


Reception of one or more signals, overlapping in frequency and time with the desired 
signal, is commonly called cochannel interference. Joint detection is the optimal minimum 
probability of error decision rule for cochannel interference. This dissertation investigates 
the optimum approach and a number of suboptimum approaches to joint detection when 
a priori information based in fields, or sets of transmitted symbols, is available. In the 
general case the solution presents itself as a time-varying estimation problem that may be 
efficiently solved with a modified Bahl, Cocke, Jelinek and Raviv (BCJR) algorithm. 

The low-entropy properties of a particular signal of interest, the Automatic Identifi- 
cation System (AIS), are presented. Prediction methods are developed for this signal to be 
used as a priori information for a joint field-based maximum a posteriori (MAP) detector. 
Advanced joint detection techniques to mitigate cochannel interference are found to have 
superior bit error rate (BER) performance than can be obtained compared to traditional 


methods. 
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EXECUTIVE SUMMARY 


Many receivers today operate in an interference-limited environment. In a dense 
environment, the performance of a receiver is limited by interference from multiple signals 
at the receiver rather than signal-to-noise ratio. In the interference-limited environment 
there is a diminishing return from investing additional resources into improving traditional 
receiver parameters such as noise figure. Advanced processing techniques exist that can 
help recover information that would otherwise be lost using a single-channel receiver [2]. 
This dissertation investigates some of these techniques. 

Interference may originate from many sources: users on the same or adjacent fre- 
quency band, unintended emissions, intermodulation, etc. This work investigates inter- 
ference from cochannel information bearing signals. Cochannel interference is commonly 
defined as the reception of two or more signals at the receiver overlapping in frequency and 
in time. 

Although multi-user detection is a mature research field [2], little prior research has 
focused on low-entropy signals. Optimal detection uses all the available information to 
make a decision. The optimal joint detection techniques are often dismissed due to the 
complexity of the receiver. Situations exist where this complexity is justified and is the 
most cost effective way to recover a signal. The field-based maximum a posteriori (MAP) 
joint detection algorithm uses the available a priori information to aid in jointly making a 
decision of what values are sent. Conceptually it is not difficult to understand; the idea is 
to select the combination of transmitted signals that maximize the a posteriori probability 
of a transmitted field. Let A be a N x N, matrix of the transmitted symbol vectors 
A = {aiQ2...a@n,}, where N, is the number of signals and N is the length of each 
vector. The optimum receiver (minimizing the probability of incorrect A) is defined as 
the MAP receiver [3], selecting the most probable A given the received vector r. Finding 


A 


the most likely set of transmitted signals, A, is now a combinatorics problem; try all 
2NNs permutations of A and select the one that maximizes the a posteriori probability. 
This method is impractical for all but the shortest of messages. This research investigates 
efficient methods of incorporating the information available at a receiver in order to make 
a decision. 

The Automatic Identification System (AIS) is a ship- and land-based tracking and 


communications system operating in the very high frequency (VHF) maritime band. The 
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primary function of AIS is to provide information for surveillance and the safe navigation of 
ships [1,4]. The AIS typically sends ship-based tracking messages indicating position and 
state information at intervals of 2—10 seconds . This dissertation uses the AIS signal as an 
example of a signal with the following properties: narrowband, cochannel and low-entropy. 
The AIS signal is shown to have significant redundancy of information content from message 
to message. Analysis of the AIS data suggests a 168-bit AIS packet contains approximately 
20 bits of new information. This research investigates the low-entropy property of the AIS 
signal and presents experimental results quantifying this property. 

This work develops the Joint Field-Based MAP algorithm in order to achieve the 
goals stated above. By incorporating a priori information from the signal sources, this 
detector out-performs maximum likelihood (ML) based joint detectors. The field error 
rate (FER) performance gain is highly dependent on the specific signal characteristics. 
Here, results are presented for representative signals that show gains of approximately 3 
— 18 dB over the current state-of-the-art. This gain does come at a cost. There must 
be a source of a priori information; this may be previous signal receptions or some other 
source of side information. Without knowledge of one or more of the cochannel signals 
there will be no gain in performance. There is also a processing cost. Implementing the 
Joint Field-Based MAP detector is both more computationally costly and more difficult to 
set up than that of the joint ML detection. The complexity of the Joint Field-Based MAP 
detector is exponentially dependent on both the number of uses and the field length. A 
complexity reduction is possible by recognizing that some field values occur with negligible 
probability. This technique significantly reduces the computational cost while incurring a 
small penalty in FER performance. 

This work focuses on advanced reception techniques that are of particular relevance 
to the military. Most of these techniques are applicable when there is sufficient extra 
processing capability available. These techniques are tailored toward asymmetric com- 
munications; the scenario where there are existing transmitters and an advanced receiver 
platform. Although this work uses the AIS signal as an example, it is envisioned that 
with adequate processing resources these algorithms can be used to improve the bit error 
rate (BER) performance of many other systems. 

This work presents simulation and analytical bounds demonstrating the benefit of 


using these advanced detection techniques on cochannel signals. 
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I. INTRODUCTION 


As the density of wireless transmitters increases, cochannel interference becomes 
the major limitation to the performance of a communications receiver. When a receiver is 
operating in an interference-limited environment additional techniques, beyond improving 
traditional receiver parameters, must be employed to increase performance. Interference 
rejection techniques such as adaptive equalization, the constant modulus algorithm, and 
various non-linear techniques can all be used to help address the cochannel interference 
problem by attempting to remove the undesired interferer [5]. These interference rejection 
techniques are not optimal. If multiple receivers are available, array processing techniques 
can be used to attenuate the undesired signal [6]. Although array processing can pro- 
vide great performance gains in a cochannel interference environment the requirement for 
multiple antennas is often impractical. This work focuses on single-antenna joint detection 
techniques [2]. With joint detection all of the signals are demodulated simultaneously. This 
method can be shown to satisfy an optimality criterion (minimum probability of error) and 
does not require multiple antennas. 

Some basic definitions should be given here. 

The term cochannel refers to two or more signals arriving at a receiver overlapping 
in frequency and time [7]. The signals may all be desired, or may consist of both desired 
and undesired signals. 

The term narrowband refers to signals with minimal excess bandwidth i.e., the 
bandwidth beyond the Nyquist bandwidth of ca where T;, is the symbol period [8]. It is 
possible to design a signaling waveform with a large amount of excess bandwidth, many 
times that of the Nyquist bandwidth, such that a single-channel receiver performs well 
in the presence of multiple cochannel signals. For example, code division multiple access 
(CDMA) is a technique that uses waveforms with excess bandwidth to enable simultaneous 
users. Other multiple access techniques, such as the time division multiple access (TDMA) 
method, typically do not use a waveform with excess bandwidth. 

The term low-entropy describes an information source that has other than a uni- 
form distribution (some symbols or sets of symbols are more likely than others) [9]. The 
low-entropy property of a signal can be exploited to achieve better bit error rate (BER) 


performance in a detector. 
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Figure 1.1: Incorporating All Knowledge into the MAP Receiver 


What is desired is an algorithm that uses the received signal and all of the other 
information available to make a decision. The problem can be formulated as follows: Let 
a, be the length N symbol column vector transmitted by the k-th user. Let A bea N x N, 
matrix of the transmitted symbol vectors A = {a1Q@2...an,}, where N, is the number 
of signals. The optimum receiver (minimizing the probability of incorrect A) is defined 
as the maximum a posteriori (MAP) receiver [3], selecting the most probable A given the 
received vector r. 


A = arg max Pr(Alr) (1.1) 
AcA 


where A is the set of possible values for A. where r is the received vector (possibly raw 
samples, possibly pre-processed by some yet to be considered demodulator). What we 


really want to maximize is 


A = arg max Pr(Alr;e) (1.2) 
AcA 


Where e is the set of other information known about the problem, as illustrated in Fig- 


ure 1.1. 


This dissertation works towards creating an efficient realization of Equation 1.3, 
ultimately resulting in 
A = arg max Pr(Alr;e), (1.3) 
AcA; 


where A; is a subset of A containing the values of A with non-negligible probability. 


A. NOTATION 


Throughout this work the following notation is used unless otherwise stated. Vec- 
tors are assumed column vectors and are represented as lower case bold face symbols, a. 
Matrices are written in upper case bold face, A. Transpose and complex conjugate trans- 
pose are written as (-)? and (-)” respectively. Statistical expectation is written as E{-}. 
An estimate of a random variable is written as er Brackets are used to index into a vec- 
tor, z[n] is the n“” element of the vector x; parentheses represent continuous time values, 
x(t). A vector with a prime, y’, indicates that it is an interleaved vector. A probability is 


represented with Pr|-]; p(-) indicates a probability density function (PDF). 


B. COCHANNEL SIGNALS 


Cochannel interference is said to occur when two or more signals are received in the 
same frequency band at the same time. Figure 1.2 illustrates a time-frequency plot with 
intensity representing the magnitude of two cochannel signals. Notice that the signals in 
Figure 1.2 overlap in both frequency and time, but portions of the signal do not overlap. 
This gives hope that the signals can be separated. 

Some communications systems such as CDMA use waveforms which are designed 
to operate when there is cochannel interference [2]. For a system using CDMA, each 
waveform is separated in code space such that the signals are orthogonal’, or approximately 
orthogonal to each other. If the signals are orthogonal then the optimum receiver may be 
realized by constructing a set of N matched filters, one for each desired signal, followed 
by N independent decision rules. In this case the optimal decision rule can be derived for 


each signal by observing the output of a single matched filter. 





‘Orthogonal is used here to mean the inner product of the underling waveforms is equal to 0, i-e., 


J diedt = 0 





Figure 1.2: Cochannel Signal Example 


Other characteristics that allow signal orthogonality (or approximate orthogonality) 
when the signals occupy the same frequency band at the same time are polarization, phase, 
and spatial diversity [2,7,8]. The focus of this dissertation is the problem of separating 
cochannel signals when the signals are not directly separable by any of the aforementioned 
methods. 


C. AUTOMATIC IDENTIFICATION SYSTEM (AIS) 


An example of a system subject to potential cochannel interference is the recep- 
tion of the ship-based AIS signal on a high-altitude platform such as an unmanned aerial 
vehicle (UAV) or satellite [10]. The AIS is a ship/land based tracking/communications 
system operating in the VHF maritime band. The primary function of AIS is to provide 
information for the safe navigation of ships and safety of life at sea (SOLAS) [1]. The 
AIS sends ship-based tracking messages indicating position at update intervals typically 
between 10 and 2 seconds. The International Maritime Organization (IMO) has ruled that 
all ships of 300 gross tons and greater must carry an AIS transponder by 2008 [1]. 

The AIS uses self organized time division multiple access (SOTDMA) to separate 
users. Each user within line-of-sight (nominally 55 km) attempts to find an unused time 
slot. Each user effectively creates his own cell centered around his ship. In general, there 
is no cochannel interference problem for users on the surface of the Earth; users of the 


same time slot are either not within line-of-sight or are far enough away for the message 





Figure 1.3: Example of Reception from Multiple Cells 


importance to be relatively low. A potential interference issue arises when the reception 
platform has a field-of-view larger than the nominal 55 km as illustrated in Figure 1.3. 
The coded cells represent areas where the SOTDMA algorithm has separated the signals 
via TDMA. Cells of the same line style indicate transmissions during the same time slot. 

The AIS is not the only signal that can benefit from joint demodulation methods. 
There are other examples of low-entropy signals that can benefit from these methods such 


as radio frequency identification (RFID), and signals from sensor networks. 


© Simultaneous Cochannel Signal Receptions 

Here, a model is developed in order to investigate the likelihood of cochannel recep- 
tion. Consider N synchronous users, all within line-of-sight of a single receiver, each with 
probability p of transmitting at time n. Let the discrete random variable C’ be the number 
of cochannel receptions at an instantaneous observation. C’ is assumed to be a binomial 


random variable with distribution [11] 


Pr(C = k) = p*(1 — pp) ({.) (1.4) 
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Figure 1.4: Probability of k Simultaneous Receptions 


This model describes a shared network of N synchronous users grouped into cells where p 
is the probability of any user transmitting at time n. 

As an example, consider 40 cells, each cell with M users, where each user has a 
probability of transmission p = 5/2250 = 0.002. This corresponds to an AIS reporting rate 
of five reports per minute. Figure 1.4 illustrates the probability of receiving k cochannel 
transmissions at any one time. The peak of each curve corresponds to the most likely num- 
ber of receptions for a particular source density. As the number of users per cell increases, 
the probability of simultaneously receiving more than one signal at a time increases. Notice 
that for a low density (low number of users per cell) the probability of cochannel reception 
is low. As the density of users goes up, the most likely number of simultaneous receptions, 
k, increases. For this example, with 37.5 users per cell, the most likely outcome is simul- 
taneous reception of three signals. The important thing to understand from this example 
is situations exist where a receiver is likely to receive a small number (less than three or 
four) cochannel receptions. Although the work in the dissertation generalizes to N users 


the focus is on two or three users. 


D. COCHANNEL SIGNAL RECOVERY TECHNIQUES 


In order to successfully jointly demodulate cochannel signals there must be some 
distinguishing feature between them. In each case the distinguishing feature leads to 
underlying signaling waveforms which have a normalized cross-correlation coefficient [2] 


with magnitude less than one. 


1. Frequency Separation From Doppler 

If the transmitters of the signals s;(t) have different relative velocities with reference 
to the collection platform they will appear shifted in frequency with respect to each other at 
the receiver. This frequency separation can be exploited in the receiver. For narrowband, 
high frequency signals with a large difference in relative velocities, it is quite possible that 
the relative motion would completely separate the signals at the receiver such that they 
are orthogonal. As an example, the space shuttle at low earth orbit (LEO) receiving the 
AIS signal would see a maximum frequency separation of approximately 4 kHz. While this 
is enough to assist in joint demodulation, it is not enough to separate the signals such that 


they are orthogonal, so single-channel demodulation techniques are not effective. 


2. Multiple Receivers 

Multiple receivers provide an additional level of diversity that may be exploited. 
Much work has been published in the area of beam-forming with multiple coherent re- 
ceivers; see [6,8] for a thorough overview. Although this is a promising area, multi- 
ple receivers result in additional cost are sometimes impractical, therefore this disser- 
tation will focus on single-antenna architectures. Single-antenna receiver architectures 
for cochannel interference mitigation are often referred to as Single-Antenna Interference 
Cancellation (SAIC) [7]. 


3. A priort Symbol Probabilities 

Another promising area is taking advantage of a priori symbol probabilities when 
performing the demodulation. If data contained in the signal s;(t) are not random but 
known or partially known, a joint receiver structure can be designed that incorporates this 


knowledge. 


Table 1.1: Binary Representation of Position 





Binary Position (m) 
00000000 +0 meters 
00000001 +1 meters 
00000010 +2 meters 
11111111  —1 meters 
11111110 —2 meters 














Consider a transmitted position value (x). The position value, in meters, is mapped 
to an 8 bit two’s complement representation (Table 1.1). The bits are mapped to binary 
symbols a € {—1,1}, and transmitted through an additive white Gaussian noise (AWGN) 
channel. This position value is transmitted via a binary antipodal modulation scheme. 
The transmitted signals are then received along with AWGN, and are coherently detected. 

First, suppose there is no a priort information available for the position. The 


probability of receiving an incorrect position P, is given by [8] 
Ie =i be (1 = tarred ee (1.5) 


where Pyiterror is given by [8] 


E 
Pie esion = Q ( #| ; (1.6) 
Q(x) is defined as 
Q(x) = os dr, (1.7) 


and Ey, is the squared Euclidean distance between transmitted symbols for binary 1 and 
binary 0 
Ea = (ao = a1)" (ag = a) (1.8) 


where ap and a, are the symbols used to represent a binary 0 and a binary 1 respectively, 


and Np is the one-sided noise power spectral density. 


Now suppose it is known that the actual position value is either 15 or 16 meters, 
x € {15,16}, with equal likelihood. The binary representation of the message therefore is 
either 
00001111 


or 
00010000 


Although these two values differ in information by one bit, the Hamming distance? in the 
“code words” is 5 bits. Not only has the problem been reduced from selecting one of 2° 
sequences to that of selecting one of 2! sequences, the sequences are further separated in 


Euclidean space, making the decision less error prone. Here, the probability of position 


n= 0( (He). (19) 


where Eq is defined as the squared Euclidean distance between the symbols representing 


error is only 


a binary 1 and a binary 0 as in (1.8). This is a significant improvement in probability of 
error over the first case. 
Now suppose the position was known to be either 16 or 17 meters, x € {16,17}, 


with equal likelihood. The positions are represented as 
00010000 


and 
00010001 


respectively. In this case, the Hamming distance between words is only one yet, by using 
this a priori knowledge of the transmitted values the decision is selecting one of 2!, instead 


of selecting one of 2° possibilities. The probability of position error is then given by 


Ea 
P=] ( ms (1.10) 


where Eq is given by (1.8). 





?Given two strings of equal length, the Hamming distance is the number of positions that the strings 
differ [12] 


There are two challenges to exploiting a priori probabilities in a receiver. The first 
is the calculation of the a priort probabilities, and then mapping this information to the 
received symbols. Estimates for physical quantities such as position, velocity, time, etc. 
can be performed with a variety of predictors. The second is the increased complexity of 
the receiver. For the example above where the transmitted value is one of two equally 
likely position values, there is little increase in complexity. For a large number of possible 


position values that are not equally likely, the decision rules become complex. 


4. Time-of- Arrival Offset 

Many communications systems such as AIS are packet-oriented. For some collection 
scenarios even packets in a synchronized network could arrive with a time offset. Over- 
lapping packets should be treated as a special case. Three situations are illustrated in 
Figure 1.5: clean reception of a packet, clean reception of the beginning or end of a packet, 
and overlap of two packets. 

Clean reception can be used for calculation of symbol probabilities. The clean 
beginning or end of a packet can be used for phase and frequency estimation. The region 


of overlap is where the joint detection techniques are used to recover the signals. 


5. Multiple Message Integration 
If a data source transmits the same information in multiple messages, soft deci- 
sions from each message can be combined to obtain a more reliable estimate of what was 


transmitted. 


E. SIMPLE MITIGATION EXAMPLE 


It is instructive introduce a simple example of cochannel interference mitigation. 


Consider two synchronous symbol streams: 


n = GO: do QF 28° 4 
a In] = {-1, +1, +1, -1, —-1} 
a) [In] = {-1, -1, +1, +1, -1} 
r[n] =a +a® {-2, 0, +2, 0, —2} 
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Figure 1.5: Packet Overlap Example 


Given the received stream r a reasonable estimate of a® and a) is 
ox) = qe) = {-1, +1, 41, +1,-1}. (1.11) 


There are two locations for possible errors in each of the received streams, the locations 
at n = 0 and n = 3. At both of these locations, if the bit for one stream (a or a@)) is 
known the bit for the other could be determined. 

Now suppose that the data sequences a and a) are sampled at twice the symbol 


rate such that: 


n = 0 1 2 345 6 7 8 9 
an =aM[[2]] = {-1, -1, 1, 1, 1,1, -1, -1, -1, -} 
an] =a [[2]] = {-1, -1, -1, -1, 1,1, 1, 1, -1, -1O 

r =a’ + a’) {-2, -2, 0, 0, 2, 2, 0, 0, —2, —2} 


This provides no more information than the first case. 
Now suppose that the symbol timing from s; and 5s is received such that s is 
delayed by T,/2 and this delay is known at the receiver. This delay could be caused by 


unsynchronized transmitters or path length differences in synchronized transmitters. 


GE 


n= 0 1 2 38 45 6 F 8 9 46 
(Pi tS, Wie oe) eae 
(0s S30c 0, es, as. Sie; a, he en SS 
De TOA 10k. Qe SO - =e 


aM In {2 


ore) = 5) 


r’=-aY+a"2) = {-1, -2, 0, 0, 0, 


a wis 
== 
l 


The two transmitted sequences may now be resolved from the received sequence. An 


algorithm for recovering a and a) is as follows. 


initialization: 
a 10] = r”[0] (1.12) 
60] = r[1] — a [0] (1.13) 
recursion: 
6 [k] = r"[2k] — dO [k — 1] (1.14) 
6k] = [2k +1) — a Ik] (1.15) 


This example can be extended to continuous time digital communications signals 
received by a bank of matched filters sampled at the symbol rate, followed by a joint 
detection algorithm. Because of the recursion, a single error at step n can lead to multiple 
errors for later steps. The propagation of errors is investigated further in Chapter V and 
Chapter VII. 


F. IMPORTANCE OF THIS RESEARCH 


Multi-user detection is a mature research field, yet there are still some areas where 
new contributions can be made; little work has focused on low-entropy signals. This work 
focuses on the reception of low-entropy signals from simple (no advanced coding techniques) 
transmitters with an advanced single-antenna receiver. The first assumption is that the 
transmitter is not generating signals intended for joint reception, and the transmitted 
signals have some degree of predictability. The second assumption is that the receiver has 


processing resources available to implement advanced algorithms. This research extends 
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prior techniques of joint reception to this scenario, allowing for successful detection of 


signals with a lower error rate than was previously possible. 


G. 


CONTRIBUTIONS 


The major original contributions of this work are the following: 


A maximum likelihood sequence detection (MLSD) for a bit stuffed data source. 
This detector gives superior BER performance compared to the ML symbol-by- 
symbol detector for the same signal-to-noise ratio (SNR) (a gain of 0.2 — 1.0 dB). 


An AIS communications state predictor and demonstration of the performance for a 
representative data set. The communications state field in the AIS message contains 
important synchronization information that can be used to predict when it is likely 


another similar message will be sent. 


Prediction method for AIS field data; including the latitude and longitude fields, 
and simulation of the performance of the prediction methods. The 27 bit longitude 
field is shown to contain approximately four bits of new information once a single 


message has been received. 


A Joint Field-Based MAP detector for cochannel signals. This detector results when 


a priori information is available for fields of data. 


A method of state reduction; by ignoring states with negligible probability (similar 
to the algorithms presented in [13,14]). State reduction becomes necessary for the 


general solution to make the problem computationally feasible. 


An analysis of time-of-arrival variation for AIS packets. Accurate estimation of 


timing parameters is an important component of the optimum detector. 


An analysis of the reporting accuracies for the state information fields in an AIS 
message, for example, the reporting precision for the longitude and latitude fields 


are found to vary from AIS to AIS. 
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H. DISSERTATION ORGANIZATION 


The remainder of this dissertation is organized as follows: Chapter II provides an 
introduction to some fundamental concepts of joint detection and a full literature review. 
Chapter III provides a complete overview of the AIS and presents some experimental re- 
sults demonstrating the “low-entropy” property of this signal. Chapter IV provides an 
overview of prediction methods and develops predictors for the AIS signal specifically. It 
provides experimental results showing the performance of the predictors on collected data. 
Chapter V develops the MLSD for a bit-stuffed data source; bounds on BER performance 
are developed and compared to simulation. Chapter VI is an overview of the multicat- 
egory classification problem that results from a field-based MAP detector based on field 
a priori information. The single-channel detector is then extended to the special case of 
synchronous cochannel reception. Bounds on performance for the single-channel and syn- 
chronous cochannel example are developed and compared with simulation. Chapter VII 
extends the results of Chapter VI to the general case of asynchronous signals. In this ex- 
tension, a time-varying trellis is developed to carry forward a priori information from each 
field. An efficient trellis traversal method is presented. Finally, Chapter VIII summarizes 
the results of this work, provides guidance for the practical implementation of some of this 


ideas from this work, and discusses areas of future research. 
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II. PRIOR WORK 


This chapter provides a history of prior work in the area of cochannel interfer- 
ence mitigation and multi-user detection. A review of other authors’ interest in AIS is 
presented. An overview of prior research in both coherent and non-coherent Gaussian min- 
imum shift keying (GMSK) receivers is developed. Finally, a review of the fundamentals 


of communications theory required later in this work is provided. 


A. AIS 


In this dissertation, AIS is the exemplar low-entropy signal considered for joint de- 
tection. There has been some interest in collection of AIS from high-altitude unmanned 
aerial vehicle (UAV) and from space [15-17]. In [17], the authors discuss the future of mar- 
itime surveillance using micro-satellites; they investigate collecting radar pulses and/or AIS 
signals along with satellite-based synthetic aperture radar (SAR) imagery. The authors also 
discuss the Norwegian micro-satellite NSAT-1 which is aimed at demonstrating detection 
and direction finding on maritime X-band radars. 

As more ships are equipped with AIS and as the next class of low power AIS (Class 
B systems) become available, cochannel interference is expected to become more of a 
problem [18]. In [17], the authors discuss the potential interference issues associated with 
collection of AIS from space. A Norwegian student project, NCUBE, is looking into the 
feasibility of collecting AIS from space. In [17], the authors hint at “advanced on-board 
processing” to enable a satellite to handle more messages than would otherwise be possible 
from a single receiver. 

In [16], the authors discuss some modifications to AIS to support Long Range Iden- 


tification and Tracking (LRIT) in an effort to increase global maritime domain awareness. 


B. GMSK 
GMSK is the modulation scheme used in the AIS. GMSK is a popular modulation 


technique when spectral efficiency and constant modulus is desired. A constant modulus 


signal requires a less complicated amplifier stage at the transmitter making this an attrac- 
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tive scheme. Murota [19], investigates GMSK modulation for what has become today’s 
Global System for Mobile Communications (GSM) cell phone systems. Murota examines 
spectrum manipulation of minimum-shift keying (MSK) by introducing a pre-modulation 
Gaussian low-pass filter (LPF). Murota explores the general properties of GMSK, error 
performance, and bandwidth consumption, and investigates some implementations. 
GMSK is a continuous-phase modulation (CPM) signal with a Gaussian filter, g(t). 


The complex baseband representation of a CPM signal is given by: 


2(t) = Aexp(je(t, @)), (2.1) 


where A is the signal amplitude, assumed to be one hereafter, and ¢(t;@) is the phase 


term whose current value depends on all previous symbols (a,), 
o(t; a) = hr S- ang(t — nT;). (2.2) 


Here, 7, is the symbol period and h is the modulation index, assumed to be a rational 
number. The modulation index h is proportional to the maximum phase change per symbol 
period. The maximum phase change per symbol period is ha/2 radians. The information 
bits are mapped to symbols a, 

a, € {-1,+1}, (2.3) 


and q(t) is the integrated pulse shape 


For GMSK g(t) is the pulse response of a Gaussian filter, 


1 = 
g(t) = rect(t/T;) ok ang? T2 exp (sa ; (2.5) 


where o = ,/In(2)/2aBT, [20], * denotes convolution, and rect(t/T;) is defined as 


1/T, |t|<T,/2 


rect(t/T,) = 
0 [t] 72 
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The pulse response g(t) can be written as a difference of Q functions!, 


t—T,/2 t+T,/2 
Q Caw —-Q (see) : (2.8) 


ae! 
OF. 





g(t) 





where BT, is the time-bandwidth product of the pulse response g(t), and B is the 3 dB 
bandwidth of the Gaussian filter. 

Typically g(t) will be truncated to LZ symbol periods, depending on the value of 
BT; such that g(t) is nonzero on the interval [0, LT] and symmetric around LT,/2. The 
integrated phase function, q(t), is nonzero on the interval |[0,00) and is constant 1/2 for 
values of t > LT, [8]. 

CPM signals are in general non-linear; they can not, in general, be represented as 
a superposition of amplitude modulated pulses. Laurent [21] shows that for finite pulse 
response, g(t), any binary CPM signal can be represented as a superposition of 24~! am- 
plitude modulated pulses. This is an extension of the well known representation of MSK 
as offset quadrature phase shift keying (OQPSK). In [22], the authors develop a minimum 
mean-square error (MMSE) pulse representation of CPM signals, and provide a supple- 
mentary proof of the Laurent representation. In [23], the authors investigate modulator 
and demodulator simplifications based on the Laurent representation when the modulation 
index h = 1/2. 

In [24], the authors show that for small modulation index, approximately less than 
0.5, the performance gain from using additional pulses after the first pulse in the Laurent 
representation is minimal. In [25], the authors compare the performance of a MLSD using 
only the “main pulse”, referred to as the AMP-MLSD. They show the BER performance 
degradation from using only the main pulse in the MLSD is a fraction of a dB. The 


importance of the Laurent representation is that only the first few terms in the sum are 
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needed, thus reducing the complexity of a receiver. A full description of the Laurent 


representation is provided in Appendix A. 


C. SYNCHRONIZATION 


Carrier phase and frequency information are required for coherent detection of 
GMSK signals. Typically this is accomplished with an initial estimate from a known 
training sequence followed by some form of phase tracking algorithm. 

Boashash [26], provides an exhaustive overview of estimation techniques for in- 
stantaneous frequency and, by extension, phase. He covers various maximum likelihood 
estimate (MLE) techniques as given in Kay [27], and model-based methods such as mod- 
eling the phase term as a linear equation and using least squares to develop an estimate 
of frequency and phase as in [28]. In [28], Tretter shows that for high SNR the frequency 
estimate based on fitting a line to the unwrapped phase term has the same variance as the 
Cramer-Rao bound. In [29], Hansen provides an excellent tutorial on phase and frequency 
estimation and provides additional examples. 

In [30], the authors investigate a method of jointly estimating both data and timing 
in the presence of intersymbol interference (ISI). The authors model the ISI as a Markov 
chain, and model the timing variation as a random walk. 

Rice, et al. [31], discuss data-aided phase estimation for GMSK using only the first 
Laurent pulse. They show that the bias of the ML estimate in [32] is due to the edge effects 
of the estimation in that paper. Mengali and D’Andrea [32] serves as an excellent refer- 
ence for synchronization techniques of sampled signals. It provides a thorough theoretical 
background without losing sight of implementation details. The authors cover both opti- 
mal and suboptimal techniques, and provide closed form solutions for performance where 
practical. 

Harris, et al. [33], discusses timing synchronization using polyphase filters. The 
authors address the problem of re-sampling and interpolating after a matched filter by 
implementing the matched filter as a polyphase filter bank. 

In [34], the authors discuss joint carrier-phase and symbol-timing for coherent CPM 
receivers. 

In [35], the authors show that the Cramer-Rao bound (CRB) for data aided timing 


or carrier phase recovery varies with training sequence. 
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There are various methods that combine symbol and channel estimation into the 
same problem. Reader, et al. [36] introduce a Blind Maximum Likelihood Sequence 
Detection (BMLSD) where the estimation of the channel parameters is included in the 
derivation of the MLSD. In [37], the authors investigate the blind acquisition character- 
istics of Per-Survivor Processing (PSP). They show that short-term poor performance of 


PSP is often from transmitted sequences that do not allow the channel to be determined. 


D. GMSK RECEIVERS 


GMSK signals may be demodulated coherently or non-coherently [38]. Huang, et 
al. [39], provide an excellent overview of coherent and non-coherent receivers for GMSK. 
They propose a fast Fourier transform (FFT) based architecture for phase, carrier, and sym- 
bol timing recovery along with demodulation both coherently and non-coherently. They 
perform simulations for receiver performance for various estimation errors in phase, fre- 
quency, and symbol timing. In [40], the authors analyze decision feedback for one-bit 
differential GMSK detectors. In [41], the authors analyze the performance of one-bit dif- 
ferential detection of GMSK signals. In [42], the authors analyze the performance of two-bit 
differential detection. In [43], the authors develop a simplified non-coherent detector that 
uses a simplified phase pulse. 

After Laurent [21] proposed representing a CPM signal as a linear combination 
of amplitude modulated pulses, many authors investigated using this method to develop 
receivers based on all or a fraction of the Laurent pulses. In [44], Kaleh uses the Laurent 
representation of a CPM signal to develop the optimal receiver as a set of N complex 
matched filters followed by a Viterbi detector. He then demonstrates that by using a 
subset of the N complex matched filters a simplified Viterbi receiver can be developed that 
approaches the performance of the optimum receiver. He then provides an example using 
GMSK with BT = 0.25. In [45], the authors develop a four-state maximum-likelihood 
sequence estimation (MLSE) GMSK demodulator. They achieve state reduction by only 
using the first Laurent pulse and by using a whitening filter prior to the MLSE. 

Many authors have looked at differential detection of GMSK using decision feed- 
back [46,47]. There are many proposed non-coherent detection methods for GMSK sig- 
nals. In [48], the authors propose a non-coherent differential detection method based on 


maximum-likelihood sequence estimation of the transmitted phase. In [49], the authors 
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investigate the performance of decision feedback equalization and non-linear data-directed 
estimation as alternatives to MLSE in the presence of adjacent channel interference. 

In [47], the authors investigate GMSK differential detectors with decision feedback 
for cochannel interference (CCI) channels. The authors show the superior performance of 
a 2-bit differential detector compared to a 1-bit or 2-bit detector with one interferer and a 


carrier-to-interferer power ratio of 11 dB. 


FE. MLSE 


MLSE is commonly used as a detection method for signals where ISI is introduced 
either at the transmitter or in the channel. The Viterbi Algorithm is an efficient method for 
finding the shortest (or the longest) path through a layered graph or trellis. The algorithm 
has been put to use as an effective method to compute the MLSE. 

In [50], Sklar provides an excellent tutorial on using the Viterbi decoding algo- 
rithm for channel equalization. He introduces the Viterbi decoding algorithm, discusses 
computational complexity and applies it to channel equalization. 

The Forney MLSE is introduced in [51]. This receiver structure is a matched filter 
followed by a whitening filter to remove the statistical correlation of the noise from sample 
to sample. The output of the whitening filter is then applied to a sequence estimator such 
as the Viterbi decoder with branch metrics of Euclidean distance. The matched filter, 
sampler and discrete-time whitening filter are often together referred to as the whitening 
matched filter [8,52]. In [53], a new upper bound is developed that is shown to be tighter 
than the Forney bound developed in [51]. 

The Ungerboeck receiver, introduced in [54], does not require a whitening filter. 
It consists of a matched filter and sampler, followed by a sequence estimator such as the 
Viterbi decoder with a modified branch metric. The Ungerboeck receiver is developed from 
the same ML principle as the Forney receiver and will have identical performance in the 
same channel. The Ungerboeck receiver uses the symmetry of the autocorrelation function 
of the channel in development of the branch metric. This modified branch metric negates 
the requirement of the whitening filter in the development of the MLSD. 

A unification of both the Forney and Ungerboeck receivers is presented in [55]. In 


this work, each receiver is derived from the other providing further insight into each algo- 
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rithm. The authors show that the whitening operation in the Forney receiver is canceled 


by the Euclidean distance metric, leaving the modified metric of the Ungerboeck receiver. 


F. COCHANNEL INTERFERENCE 


There are many terms used in the literature that generically refer to cochannel 
interference. The term multiple access interference (MAI) is commonly used to describe 
interference from users within the same network or cell. CCI is often used to describe inter- 
ference from a different network or different cell. In this dissertation cochannel interference 
is used to describe both phenomena. 

Shnidman [56], was the first to point out that cochannel interference and ISI are es- 
sentially an identical phenomena. Both CCI and ISI lead to multiple symbol contributions 
at the output of a filter matched to a single symbol period. Kaye [57], further investigated 
the problem of multiple signal detection over multiple channels. 

Single-antenna cochannel interference mitigation techniques can roughly be placed 
into two different categories. The first category consists of filter based approaches including 
subspace methods that involve filtering out the interference. These non-linear techniques 
rely on the signal and interferer being somehow linearly separable. The second category 
consists of multi-user detection techniques such as joint maximum likelihood sequence 
estimation (JMLSE) and successive interference cancellation. These techniques jointly 
estimate all of the signals; separation is achieved through the joint estimate. The authors 


in [7,58,59] provide a good overview of these techniques. 


1. CCI Filtering Methods 

There has been much work in the area of reducing the complexity of the optimal 
joint detector. The computational complexity of the optimal detector grows exponentially 
in both ISI length and number of signals [2], making it practical only for a small number 
of simultaneous users. 

In [60-62], the authors develop an interference cancelling matched filter by using 
known training sequences. In [63], the authors investigate the cancellation of unwanted 
interfering signals based on second-order statistics. The constant-modulus algorithm and 
the phased-locked loop are extended for use in blind multi-user detection in [64]. In [65], the 


authors provide a performance comparison of coherent and limiter-discriminator GMSK 
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detectors in the presence of a single interferer for various carrier-to-interferer ratios and 
various bandwidth-time products. 

Pukkila, et al. [66,67], investigate SAIC for GSM where the authors assume the 
interference is constant modulus, and they assume they do not have channel state infor- 
mation. This technique is known as the Constant Modulus - Single-Antenna Interference 
Cancellation (CM-SAIC) algorithm. 


2. Joint Detection 

Optimum joint detection techniques are often rejected because the complexity of 
the computation is exponential in the number of users. For detection of only two or three 
users the optimal techniques become viable. For narrowband signals, SNR limitations 
would typically limit the joint detection of multiple users to no more than two or three 
users. 

In [68-70], the JMLSE is developed for continuous-phase modulated signals. This 
involves representing the CPM signals using the Laurent representation [21], and following 
the methods developed for linear modulation methods. 

Polydoros and Chugg [71] discuss joint reception of narrowband signals based on 
work in [2]. The authors investigate iterative-detection techniques and spend some effort 
discussing complexity-performance tradeoffs. 

In [72], the authors develop a lower bound for bit error probability for the JMLSE 
with perfect channel knowledge. 

Van Etten [59] describes a method of extending the Viterbi algorithm for ISI to 
handle multiple cochannel signals. He termed this the “Vector Viterbi Algorithm” . 

Verdti [73] uses the Ungerboeck formulation of the MLSE to form an optimal k-user 
multi-user detector for asynchronous signals. 

In [74], the authors describe a blind adaptive joint detector. The detector jointly 
estimates both channel and transmitted data for each signal. The authors propose a T,/2 
sampler. Their detector is based on [75]. 

Giridhar, et al., [76] describe an optimal method for jointly demodulating cochan- 
nel signals termed the joint maximum a posteriori sequence detector (JMAPSD). They 
develop an optimal joint demodulator and investigate the performance of suboptimal re- 


celvers. 
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In [77], the authors develop the optimum joint MLSE for GSM signals. To enable 
coherent detection, the authors develop a method to jointly estimate channel parameters 
based on a known training sequence. 

In [78], the authors develop a model-based method for CPM demodulation in the 
presence of cochannel interference. The authors use modern spectral estimation techniques 
on a symbol-by-symbol basis in order to estimate instantaneous frequency. They provide 


an excellent summary of standard CPM demodulation techniques. 


3. Suboptimal Joint Detection 

In [79], a receiver structure consisting of a bank of matched filters followed by a 
whitening filter and a vector Viterbi algorithm is proposed. In [80], the authors propose 
an orthogonal multi-user receiver. The proposed receiver does not require knowledge of 
signal amplitude or noise power. Ranta, et al., [81] propose a method of exploiting the 
independently fading multipath channels to provide separation of transmit waveforms to 


be used in joint demodulation. 


4. Development of the JMLSE 

Below is an introduction to JMLSE. It is well known that the output of a single 
matched filter generates a sufficient statistic for detection of a signal in Gaussian noise [8]. 
Consider the detection of multiple synchronous signals in noise with a single receiver. The 


received signal r(t) is the summation of these multiple signals plus noise: 


i= 3 si,(t) + n(t). (2.9) 


A simple strategy for detecting the multiple signals is to form a bank of matched filters 
each followed by a decision threshold. If the signals are orthogonal and the symbols equally 
likely, this structure could be shown to minimize the probability of error. If the signals are 
not orthogonal, this method will not minimize probability of error. For signals that are 
highly correlated or slightly correlated with large differences in received power, this method 
will perform very poorly. The outputs (y1, y2,...,yn,) of N, matched filters matched to 
each of the signals are sufficient statistics for estimating the symbols (a1, @2,...,@,), but 
yi is not a sufficient statistic for estimating a, [2]. This is because, in general, the signals 
S$, are not mutually orthogonal, so in general the decision statistics y;, are not independent, 


even when n(t) is AWGN. It is for this reason that joint detection is required. W. van 
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Etten [82] demonstrates that the joint detection problem is an extension of the ISI problem. 
He was the first to show that the outputs of multiple matched filters matched to each of the 
cochannel signals form sufficient statistics for joint estimation of the transmitted signals. 
For asynchronous signals, the matched filter outputs at a particular sample time 
are not sufficient statistics for determining the symbols transmitted at that same instance; 
the matched filter outputs from the entire duration of the signal (yi, y2,...,yn,) must be 


used in the estimation of the symbol vectors (@1,Q@2,...,@n,) [2]. 


5. General N, Signal Case 

Below we the joint MLSD is developed for the N,-user case using a complex base- 
band representation and a modified metric introduced by Ungerboeck in [54] for an ISI 
channel. The development closely follows the work by Verdt [2, 53,73], Murphy [68] and 
examples in [8]. 

Consider the reception of N, signals in AWGN 


Ns 
r(t) = So sp(t) + 2(t). (2.10) 
k=1 
Where z(t) is AWGN and s;(t) is represented as 


sp(t) = s Axorg[i]e(t — iT — Te), (2.11) 


where the unit-energy signal pulse, w(t), is only nonzero on the interval [0,7], 7, is the 
time offset, and A;, is the amplitude for the k” signal. 


Form a set of matched filters whose outputs at the sample time are 


yt] = / r(t)y,(t — iT, — T,) dt, (2:12) 
where Ww, are the unit-energy signal pulses from (2.11). The matched filter outputs, 
(Yo, Y1;---;YN,—1),; are a linear combination of all the signals’ individual responses to 
the filter: the desired response plus the response of 2 x (NV, — 1) other symbols from the 
other channels. For pulses with a duration of one symbol period, the matched filter outputs 
will consist of contributions from the symbol just before and the symbol just after from 


each other channel. This is shown in Figure 2.1. The output of the matched filters can be 
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written as 
Uk [7] = Apa |i] 


dis > Aja; [t] We(t)b;(t — (7) — Te) dt 


j<k 


+DAvajli+ f “Un(tuy(t —T, — (7j — T)) dt 
j<k (2.13) 
+S 0 Ajajli—1] J de(t)yj(t + Ts — (7) — Te) dt 


j>k 


af ys Aja;{i] Welt) y(t — (73 — Te)) dt 


j>k 








T Zk [z], 


where z;[i] is a complex Gaussian random variable with variance 0? = No. The Gaussian 
random process formed by interleaving the noise components from the output of each 


- 
matched filter is given by z|7] = EB [¢] zo[?] ... zy, [éi]) and has autocorrelation matrices 


oR’ [1], iff =i+1 
£{zli2"]} o?R'(0], iff =i 
ee oR], iff =i-1 ae 


0 otherwise 


where Ri — j] is defined as the matrix with elements n,m given by 
Romli-a= f eh 4a rah (2.15) 


It is convenient to group the i” terms in (2.13) to obtain a vector form of the 
equation 
y|?] = R*[1JAeali +1) + ROJA,,a[] + RAAneli — 1) + nf], (2.16) 


where A,, is a diagonal matrix containing the signal amplitudes, A,, = diag{Ai,..., Ax}. 
R{0] is the cross-correlation matrix between the j and k" signal for the current symbol, 
and Ryx[0] = J, ° ve(t),(t— (7) — Te)) dt. The matrix R[1] contains the cross-correlations 
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T3 a3(0] a3{1] a3[2] 






































Figure 2.1: Three Asynchronous Cochannel Signals 


between the 7 signal at the current symbol and the k’” signal at the next symbol for 
Vie 
0, ify 25K 
R;x(1] = - : (210) 
j{i — 1 ie oe (t)yj(t+ Ts; — (7) —T)) dé iff <k 
For equally likely symbols in AWGN the MAP receiver is equivalent to the ML 


receiver that selects the symbol matrix A such that the following equation is satisfied 


A=arg min / Ir(t) — S© sp (t, en)/? dt, (2.18) 


cas k 


where the symbol matrix A is a (N, x N) matrix with elements a;[i] € {—1,1}. Equa- 
tion (2.18) can be expanded and the terms that do not depend on a, dropped, to show 


this is the equivalent to maximizing 


2 


dt 


Q(A) = 2 f r(t) Dott Qk) dt — [. 2, sets) 


=2f r(t) sults on) dt - a DD sult on) sie(t, ow) dt 


oo py 


8 r(t) D_ sult, on) ar—2 fo TE Rlowtanjst(taw)) at 
-[ S- |si(t, erg) /? dt 


Where the last step follows from using the conjugate symmetry of the autocorrelation 








(2.19) 


matrix to write the double summation as twice the real part of the sum of the lower 
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triangular plus the diagonal terms. Equation (2.11) can be used to expand the first term 
of (2.19) and (2.12) to write this in terms of the matched filter outputs 


a r(t) S> sp(t, arm) dé = 2 f r(t) S> DS Apap [i] Vx (t — tT, — 7) dt 
= 20> Anculilanli . 


(2.20) 


The second term can be expanded using (2.11) to get 


[ DD 4 (t, Qz) sii( (t , Ap) dt 


co kl! 


Fe [ pee Ds Axorsli} u(t — T's — Tr) x Apc lie (t — iTs — Te) dt. (2.21) 
=> = = Apart] An ow [J] Rin [i — J] 


Where the matrix Ri — j] is defined as in Equation 2.15. Equation (2.20) and (2.21) can 


be combined to obtain: 
= = 2D) Anal vag lp 8p Bed [t] Apo [J Rr, x [i —j] (2.22) 
The above can be written in vector format as 


A) =2)/a"[i]Am fap 2 iJAmR{i — j]/Amalj]. (2.23) 


The symbol matrix A is a (N, x N) matrix with column vectors a,. Maximizing 
Q(A) is a combinatorial problem; find A by trying all 2*" permutations of A. Using this 
method computation is exponential in both message length N and number of users, Ng, 
making it an unattractive solution. 

The key to reducing the complexity of the solution is to recognize that R[i — J] is 


nonzero only for |i — j7| < 1 due to the finite support of 7. This allows us to write 


A=arg max DRO) (2.24) 


ZF 


where #(-) indicates the real part and 4; is given by 
Ni = 207 [i] Amy] — aT [i] AmR[0| Amoi] — 207 [i] AmR[1]Amali — 1] (2.25) 


The last term has a factor of two from the symmetry of Rj,[1] = R,j;[-1]. The alija[i+ 1] 
from the previous term is absorbed into the a|i]a|i — 1] from the current term. 

Writing the maximization problem in the form of (2.24), where \; depends only 
on @[i] and ali — 1], allows us to remove the exponential dependence on N. This max- 


imization may be found using the Viterbi algorithm with 2?%s—! 


states and N stages. If 
T 

we were to form a vector y’ = lv [0] y2[0] ... yw,[O0] yi[l] ...yws[N]} as we did in 

the previous section we could perform this maximization using the Viterbi algorithm with 


2Ns—! states and N, x N stages. 


G. MAP ALGORITHM 


There are two common varieties of the MAP algorithm. Bahl, Cocke, Jelinek and 
Raviv (BCJR) [83] is a forward-backward recursion to calculate the a posteriori probabil- 
ities of the state transitions of a Markov source observed through a discrete memoryless 
channel. The BCJR algorithm is detailed in Appendix B. Abend and Fritchman developed 
an alternative algorithm [75] that does not require a backwards pass, making it suitable 
for data that is not in finite blocks or for large blocks with a time-varying channel [84]. 
Li, et al. [84] provides a thorough discussion of the two types of MAP algorithms. In 
addition, many authors have investigated efficient hardware implementations of the BCJR 
algorithm [85]. 


1. Reduced Search Techniques 

The MAP algorithm has complexity that is exponential in the number of users and 
the ISI channel length. There has been some work on reduced search techniques to lower 
the computational complexity of the BCJR algorithm [13, 14, 86]. These techniques all 
involve ignoring working probabilities that fall below a selected threshold. The methods 
can be roughly placed into two categories. The “M-BCJR” algorithm selects the M largest 
components at each stage during each pass (forward or backward) through the trellis. The 
“T-BCJR” algorithm ignores components at each stage that are below some threshold 7’. 


These techniques are geared toward their use in Turbo decoding where the BCJR algorithm 
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is used over multiple iterations. The major reduction in computational complexity is for 
later iterations. 

Vithanage, et al., [87] develop a reduced complexity BCJR type algorithm based on 
propagation of expectation . The reduction in complexity is similar to that of the M-BCJR 
algorithm proposed in [13,14], by selecting a discrete distribution with limited support. 
They term this algorithm EP-MBCJR. 


H. CHAPTER SUMMARY 


This chapter provides required background and a survey of previous work in the area 
of cochannel interference mitigation and multi-user detection. The next chapter describes 
the AIS signal in detail and provide some analysis of the low-entropy property of this 


signal. 
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Ill. AUTOMATIC IDENTIFICATION SYSTEM 
(AIS) 


A scenario subject to potential cochannel interference is the reception of the ship- 
based AIS signal on a high-altitude platform such as an UAV or a satellite. There has 
been some interest in collection of the AIS signal from both high altitude UAVs and 
from satellites for maritime traffic monitoring [10,16]. The AIS signal is particularly 
interesting for this dissertation due to its low-entropy property; the transmitted fields of 
the signal contain similar information from transmission-to-transmission. This property 


will be exploited to aid in joint detection; leading to improved FER performance. 


A. AIS OVERVIEW 


AIS is a ship- and land-based tracking and communications system operating in the 
VHF maritime band. The primary function of AIS is to provide information for surveil- 
lance and the safe navigation of ships [1,4]. The AIS sends ship-based tracking messages 
indicating position and state information at intervals typically between ten and two sec- 
onds. The International Maritime Organization (IMO) has ruled that all passenger ships, 
cargo ships greater than 500 gross tons, and all ships greater than 300 gross tons on an 
international voyage must carry an AIS transceiver by July 2008 [88]. 

An example of the end use of an AIS transceiver is shown in Figure 3.1. Figure 3.la 
shows an illustration of the output of a ship-board AIS display to aid in safe navigation. 
Figure 3.1b illustrates an aerial view of the same scene. The display provides the captain 
of a ship with a real-time picture of the local ship traffic around the vessel. The displays 
can be integrated with radar data to provide a integrated picture of local sea traffic. 

Each AIS transmits position reports at a rate based on the current state of the ship. 
For example, a ship underway and changing course is required to report more often than 
a ship that is moored. Nominal reporting intervals for Class A and Class B equipment are 
given in Table 3.1 and Table 3.2 [1]. Class A Equipment is the mandatory equipment for 
ships required by the IMO to carry AIS equipment. Class B equipment is for ships not 
required by the IMO to carry equipment. 


dl 





(a) AIS Display (b) Aerial View 


Figure 3.1: AIS Display Example 


Table 3.1: Class A Ship Reporting Interval 





Class A Ships Information / Condition Reporting Interval 








Static Information 360 sec 
Voyage Related Information 360 sec 
Safety Message as required 
Ship at Anchor moving less than 3 knots 180 sec 
Ship < 14 knots not changing course 10 sec 
Ship 14 — 23 knots 6 sec 
Ship > 14 knots changing course 2 sec 
Ship > 23 knots 2 sec 





Table 3.2: Class B Ship Reporting Interval 





Class B Ships Information / Condition Reporting Interval 








Ship moving less than 2 knots 180 sec 
Ship 2 — 14 knots 30 sec 
Ship 14 — 23 knots 15 sec 
Ship > 23 knots 5 sec 
Search and Rescue (airborne) 10 sec 
Aids to Navigation 180 sec 
AIS base station 10 sec 
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Table 3.3: AIS Transmitter Parameters 














Transmit Parameter Low Setting High Setting 
Power (W) 2 + 20% 125 20% 
Time-bandwidth product BT, 0.3 0.4 

Modulation index h 0.25 0.5 

Bit rate 9600 + 50 ppm 9600 + 50 ppm 
AIS 1 Frequency 161.975 + 3 ppm 161.975 + 3 ppm 
AIS 2 Frequency 162.025 + 3 ppm 162.025 + 3 ppm 
1. AIS Transmitter Characteristics 


An AIS transmitter on the open seas will transmit a GMSK signal with a modulation 
index h = 0.5 and transmit power of 12.5 Watts. In terrestrial waters these settings may 
be changed by the appropriate authority [1]. Table 3.3 lists the various AIS transmitter 
configurable parameters [1]. On the open seas the “high” settings are used. In terrestrial 
waters an authority can instruct the transceiver to change any or all of the parameters to 
the “low” setting. The transceiver will default to alternate transmissions between the two 
AIS channels, AIS 1 and AIS 2. 


2. AIS Multiple Access Method 

The AIS is a multi-user system. The system allocates 2250 transmission time slots 
per minute on two separate frequencies allowing for a total of 4500 transmissions per 
minute. The AIS operating in continuous and autonomous mode uses SOTDMA to sepa- 
rate users. Each user within line-of-sight (nominally 55 km radius) will monitor AIS traffic 
for a full frame (1 minute) to create a list of candidate time slots for transmission. The 
user will then select randomly among these candidate time slots for transmission during 
the next frame. The same time slot will be used for up to seven subsequent transmissions. 
Each user effectively creates his own cell centered around his ship. 

A minimal cochannel interference problem exists for users on the surface of the 
Earth; in general, users of the same time slot are not within line-of-sight. The time slot 
selection algorithm is designed such that the probability of two ships selecting the same time 
slot is low. A potential interference issue arises when receiving the signals on a platform 
that has a field-of-view larger than the nominal 110 km as illustrated in Figure 3.2. The 
three circles in the figure represent areas where the SOTDMA algorithm has separated the 
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Figure 3.2: Example of Reception from Multiple Cells 


signals via separate time slots. Within each region there is no contention; there is no more 
than one transmission per time slot per channel. At a remote receiver, like the satellite 
depicted in Figure 3.2, there is the possibility that some of the time slots will be occupied 
by multiple users leading to cochannel interference at the receiver. 

The AIS does not directly address the hidden node problem. As an example of the 
hidden node problem, consider three ships Ship A, Ship B and Ship C all on a line and, 
separated by half the nominal cell size, Figure 3.3. Ship B can see both Ship A and Ship C 
but Ship A and Ship C’ can not see each other. If Ship A and Ship C’ choose the same 
time slot, then each will interfere with Ship B. This case does not cause a problem in 
practice. If the distance between ships is large, as is the required condition for the above 
interference, the importance of receiving the message is low. 

The AIS specifies 2250 time slots each minute for each channel. Each slot is 26.67 
ms (256 bits) long. The structure of the AIS message format is shown in Figure 3.4. This 
message structure is based on the High-level Data Link Control (HDLC) protocol specified 
in [89]. The “ramp up” period allows for the transmitter to get up to full power. The 24 bit 
training sequence is a sequence of alternating ones and zeros used for synchronization at 
the receiver. The start flag and the end flag are the HDLC flags indicating the start and 
end respectively of the 168-bit data packet. The frame check sequence (FCS) is a 16 bit 
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Figure 3.3: Hidden Node Example 
































Ramp up | Training sequence | Start flag | data FCS End flag | Buffer 
8 bits 24 bits 8 bits 168 + stuff bits | 16 bits | 8 bits 24 bits 
26.67ms 


Figure 3.4: AIS Packet 


cyclic redundancy check (CRC) on the data portion of the message to increase confidence 
that a message is received without error. The CRC uses the polynomial specified in [89] 
with bits preset to one. The buffer allows for timing offsets from bit-stuffing, propagation 
delay, repeater delay, and synchronization jitter. 

Table 3.4 is a paraphrased description from the International Telecommunication 
Union (ITU) recommendation [1] of the contents of the data portion of the AIS packet 
described in Figure 3.4. There are many varieties of AIS messages. Table 3.4 illustrates 
some of the more common messages. Careful inspection of Table 3.4 reveals that many of 
the fields are either constant or predictable. The ability to predict field values is the key 


feature exploited in this work; this is explained in more detail in Section III.D. 


B. EXAMPLE AIS LINK BUDGET AT LEO 


It is useful to calculate the anticipated SNR at the receiver for an example collection. 
If the reception is severely noise-limited for single-channel reception, little will be gained 
from investigating methods to mitigate cochannel interference. An example AIS link budget 
calculation for a collector with a path length of 800 km is shown in Table 3.5. This 


calculation does not take into account factors such as multipath and considers a rather 
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Table 3.4: Type 1, 2 or 3 AIS data message break out after [1] 
























































Parameter Number | Description 
of bits 

Message ID 6 1, 2, or 3 for this message 

Repeat Indicator 2 Indicates number of times message has 
been repeated 

User ID 30 MMSI?* number 

Navigation Status 4 0 under way, 1 at anchor, 2 not under 
command ... 

Rate of Turn 8 € [—127, 127], -128 indicates not avail- 
able 

SOG 10 Speed over ground € |0, 1022] 1023 not 
available 

Position accuracy 1 1 = high (< 10m) 

Longitude 28 0x6791ACO indicates not available 

Latitude at 0x2412140 indicates not available 

COG 2 Course over ground, tenths of a degree 
€ [0,3600], 0xE10 indicates not avail- 
able 

True heading 9 Degrees € [0,359], 511 indicates not 
available 

Time stamp 6 UTC? seconds € [0,59], codes € [60, 63] 

Reserved 4 Reserved 

Spare if should be set to 0 

RAIM-flag 1 RAIM (Receiver autonomous integrity 
monitoring), 1= RAIM in use 

Communications state | 19 Specifies slot position and other TDMA 
parameters 

Total 168 




















*Maritime Mobile Service Identity (MMSI) 
universal time (UTC) 


simple gain for the transmit and receive antennas. The path loss is calculated as 


2 
Path Loss = (=) 
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where R is the range and \ is the carrier wavelength given by 3 x 10°(m/s)/162(MHz). 


The received power is given by 
Received Power = P:G;G,/Path Loss (32) 


where P; is the transmitted power and G; and G, are the transmit and receive antenna 
gains respectively. The energy per bit (£») is given by the Received Power x T, where T; is 
the bit period. The one sided noise density No is given by No = kT where k is Boltzman’s 
constant and TJ’ is the noise temperature. 

Table 3.5 demonstrates that for nominal reception from a high altitude receiver, 
the signal is not noise-limited, therefore there may be opportunities for joint reception 


techniques. 
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Table 3.5: Example AIS Link Budget 





























Parameter Val dB 
Path Length R 800.00 (km) 
Transmit Frequency f 162.00 (MHz) 
Path Loss Ls 134.69 (dB) 
Path Loss Dg 134.69 (dB) 
Power Transmitted P, 12.50 (W) 10.97 (dBW) 
Transmit Antenna Gain Gt 1.00 0.00 (dB) 
Receive Antenna Gain G, 1.00 0.00 (dB) 
Received Power ie 4.24—18 (W) -123.72 (dBW) 
Received Power P,. 4.24-13 (W) -123.72 (dBW) 
Bit Period Ty 0.10 (mSec) —-39.82 (dB-sec) 
Energy per Bit Ey AAR TS) -163.55 (dBJ) 
Representative Noise Temperature = T’ 2910.00 Kelvin 34.64 (dBK) 
Boltzman Constant k 1.3872 JK 12998 60 (dBIK*) 
Two Sided Noise Density No -193.96 (dBJ) 
Energy per Bit Ey -163.55 (dBJ) 
Two Sided Noise Density No -193.96 (dBJ) 
Energy per Bit / Noise Density £F,/No 1100.22 30.41 (dB) 
Received Power ie 4.24-13 (W) -123.72 (dBW) 
Noise Density No -193.96 (dBJ) 
Bandwidth 8 25.00 (kHz) 43.98 (dB) 
Signal to Noise Ratio SNR = 422.48 26.26 (dB) 
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C. SINGLE-CHANNEL DEMODULATION 


Before delving too far into demodulation of cochannel signals, it is instructive to 
consider a simple case of non-coherently demodulating an AIS packet in the presence of 
AWGN. This example, although not optimal in any sense, illustrates the characteristics of 


the signal. The test signal has the following characteristics: 
e data rate of 9600 baud 
e center frequency of + 160 MHz (AIS1 161.975 MHz, AIS2 162.025 MHz) 
e GMSK modulated with BT, = 0.4, h = .5 
e Frames on UTC minute boundary with 2250 time slots 
e Each slot 26.6 ms 


e Before modulation, the data are encoded non-return-to-zero inverted (NRZI) (level 


changes when a ‘0’ is sent, level stays the same when a ‘1’ is sent)". 


e The data field and FCS are bit-stuffed (when a string of five consecutive 1’s is 


encountered a ‘0’ is inserted). 


e A 24 bit training sequence of alternating 0’s and 1’s. It is not specified if the pattern 
starts with 0’s or 1’s since the data will be NRZI encoded. 


e The data contains the HDLC start and stop flags 0x7E, where the 0x notation 


indicates a hexadecimal value. 


Assume the signal has been received through a band-pass filter and mixed down 
in quadrature so the signal is analytic and centered at some intermediate frequency (IF). 
The IF is then sampled at some rate greater than twice the signal bandwidth. The first 
step in demodulation is to identify the start of the packet in time. An energy detection 
routine can be used to do this. Figure 3.5 is a time frequency plot of a single AIS packet 


in noise located by using an energy detection method. 





'The definition of NRZI is found in many sources defined as a level change when a ‘1’ is sent and no 
change on a ‘0’. The definition above for the AIS is from the ITU recommendation [1]. 
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Figure 3.5: Single AIS Signal in AWGN 


The next step is to pass the signal through a band-pass filter, removing out-of-band 
noise, and then shift the result down to baseband. Let the received analytic sequence be 
represented by 

r[n] = s[n] + 2[n]. (3.8) 


The result from filtering and mixing down to baseband is 
yln] = (r[n] * hfn]) er", (3.4) 


where w, is the approximate intermediate frequency and h[n] is the impulse response of 
the filter. The differential phase term, Ad, is given by 





Ag[n] = ¢[n] — ¢[n — I], (3.5) 
where (yln) 
Ss tage t | AU 
do) = a0 Sa mae 


and tan~'(-) is the four-quadrant arctan function. Figure 3.6 shows A@[n] for an example 
received AIS signal. The training sequence of alternating ones and zeros can clearly be 
seen. Notice that the alternating ones and zeros show up as alternating pairs of ones and 


pairs of zeros from the NRZI encoding. For example the sequence “010101010101” will be 
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Figure 3.6: Differential Phase 


either “001100110011” or “110011001100” after NRZI encoding, depending on the initial 
state of the encoder. 

The AIS signal has a training sequence of alternating 1’s and 0’s followed by the 
HDLC start flag of Ox7E. The combination of these sequences can be used to locate the 
start of the packet and to derive bit timing. 

A data sequence is then generated by decimating A@[n] with an offset and making 
hard decisions. 

Anrzi[n] = u(Ag[Dn + k] — Ad) (3.7) 


where D is the down-sampling factor, k is an integer offset, Ad@ is the mean of Ad, and 
u(-) is the unit step function. 
1 ¢>0 
u(¢) = (3.8) 
Os 62.0 
Note that Equation 3.7 assumes down-sampling by an integer. If the end goal is to recover 


symbols at some rate R,, the sample rate should be fixed at F, = 1R, for some positive 
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integer /. If this were not the case two options exist: re-sample the original data such that 
F! = 1R, or select the sample which is closest to the ideal sample time. 
Since the data ay,;z;|n] is encoded NRZI it must be converted to non-return-to- 
zero (NRZ) defined by: 
a,[n] = Gnrzi[n] ® Gnrzi[n — 1], (3.9) 


where © denotes the logical XNOR operation. Before the data are transmitted it is bit- 
stuffed. The bit-stuffing operation inserts a zero after a consecutive string of k ones. This 
has a benefit of aiding timing recovery and removing the ambiguity with detection of 


control flags. The data recovery operation is written as an operator Bz’, for example 
a = B;'(a,). (3.10) 


Where a, is the bit-stuffed vector, a is the recovered original data vector prior to bit- 
stuffing. The operator Bz! indicates the unstuffing operation when k = 5, as is the case 
for AIS. 

Data are sent byte-wise with the least significant bit first. The data are then in 
a format ready to be parsed to recover the internal data fields. A full description of all 


packet formats is given in [1]. 


D. PREDICTION OF AIS PARAMETERS 


Suppose f is a discrete random variable representing a number coded with k bits. 
If the probability mass function (PMF) of f is known to the receiver, the receiver may use 
this information in deciding which of the possible values, {0,1,...,2*~'}, of f were sent. 
When the PMF is not uniform on [0,2*~*] there will be some gain in performance by using 
MAP detection vs. ML detection. The details of the gain in performance are provided in 
Chapter VI and Chapter VII. The predictability of various discrete random variables in 
AIS data fields are discussed below. 


1. Prediction of an AIS Packet Time-of- Arrival 
Time-of-arrival prediction is broken into two estimates. The first is an estimate of 
which time slot a particular AIS transceiver will transmit on. This estimate is based on the 


SOTDMA time slot selection algorithm. The second is an estimate of the offset from the 
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Figure 3.7: Histogram of Time-of-Arrival Variation for Several AIS Transceivers 


nominal start of that time slot. This estimate is based on the nominal path length to the 
receiver and an estimation of the stability of the transmissions from the AIS transceivers. 

A particular AIS transmitter will reserve a fixed number of time slots ina SOTDMA 
frame depending on the desired transmit rate. Typically between 6 and 30 time slots per 
frame will be reserved, corresponding to update rates of ten and two seconds respectively. 
It is important to predict with some degree of accuracy on which time slot a particular 
user will transmit. This is particularly important for identifying the start of a packet in a 
cochannel reception. 

The same time slot will be used for up to seven frames before the next assignment. 
The number of frames for which a time slot is reused is determined by the SOTDMA 
algorithm. Chapter IV describes the time slot assignment algorithm in great detail. 

The specification in the ITU recommendation [1] allows for +3 symbols of jitter for 
the packet start time relative to the nominal start time. Experimental results have shown 
that under typical conditions this jitter interval is on the order of 3-100 jus; this corresponds 
to 3-25% of a symbol period. Figure 3.7 shows the distribution of time-of-arrival offsets 


for several transceivers. 
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2: Prediction of AIS Packet Data 

We may ask, “How much information is transmitted by an AIS message?” To 
answer this question the information source is modeled as a stochastic random process. 
The entropy of an error term can then be calculated, where the error term is the difference 
between the measured value and the value modeled by a random process. Most of the 
fields within the data portion of an AIS message are very predictable. It is instructive to 
analyze each field of one of the message types and present results for each. 

Analysis of data received by the United States Coast Guard (USCG) shows that 
type 1, 2 and 3 messages make up approximately 98% of AIS message traffic. Analysis 
was performed by obtaining 24 hours of coastal data from the USCG-provided application 
AlSUser and counting the number of each message type. With this in mind the analysis in 
this work is restricted to these message types; other message types are treated as unknown 
interferers. 

Four types of predictors are used depending on the field. The latitude and longitude 
fields are predicted with a Kalman filter. The communications state field uses a special 
predictor based on the protocol for SOTDMA. The time field uses an external time source 
as a predictor. All of the other fields use a first order predictor. A full description of each 
predictor along with the parameters used for each field is provided in Chapter IV. 

Some of the fields vary depending on the particular AIS transceiver used, and this 
warrants additional explanation. For example, the latitude and longitude fields are 27 and 
28 bit respectively with the least significant bit corresponding to 1/10,000 minute. Some 
AIS report to 1/1,000 minute or even 1/100 minute, however this is taken into account by 
the predictor developed in the next chapter. 

Each field is modeled as a stochastic process f(t). For each field past measurements 


at time t — 7; are used to estimate the current value at time t. 


A 


f(t) = 9(f(t—7), f(t —72),..., f(t — TN), T1, T2,---, TN) (3.11) 


The random process f(t) is discrete; the possible values of f(t) are drawn from a finite set. 
In order to measure the entropy of the probability distribution of the error in the 
prediction of f(t), the function e(t) is defined as the difference between the estimated value 


and the actual value. 
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e(t) = f(t) — fl) (3.12) 


Entropy is a natural measure of uncertainty in the random process. The entropy 
of the random process f(t) could be considered, but it is clearly non-stationary. Although 
e(t) lacks true stationarity, f(t) could not even be approximated as a stationary process. 

The entropy (in bits) of a discrete random variable p with probability mass function 
(Pi, P2,---,;PK) is given by: ' 

H(p) = S— pe logs(1/pr), (3.13) 
k=1 
where we define 0 log,(1/0) = 0. 

To find the entropy of e(t) the PMF of e(t) is required. The PMF is estimated 
using a normalized histogram, implicitly assuming that the process e(t) is ergodic. With 
the PMF of e(t), Equation (3.13) can be used to calculate the entropy (information content) 
of the field. 

Figure 3.8 shows the number of bits in each field of an AIS type 1, 2, or 3 message, 
along with the calculated entropy of the residual error process e(t). The message data used 


to create this graph is from the USCG collection platforms along the US coastline. 


E. CHAPTER SUMMARY 


This chapter introduces the AIS signal as a good candidate for low-entropy cochan- 
nel mitigation techniques. The characteristics of the transmitter are introduced and mul- 
tiple access scheme SOTDMA is discussed. The case for cochannel interference as seen by 
a remote observer is presented. Much detail is provided on the format of the transmitted 
messages and the data contained within each message. Finally measured results for the 
estimated entropy of the transmitted data is presented, illustrating the AIS signal is low- 
entropy. The next chapter discusses the specific methods used for prediction of the AIS 


messages for this research. 
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IV. PREDICTION 


This dissertation develops an algorithm for joint detection of two or more signals 
exploiting a priort knowledge of the data source for one or more of the signals. The 
fundamental assumption is that we are able to predict in some way future values of the 
desired signals, and use these predicted values to aid in demodulation. A distribution is 
desired for all possible values of a signal based on past (prediction), and perhaps both past 
and future (smoothing) measurements. This chapter introduces some methods of prediction 
and applies them to the AIS signal. Results are then presented for the application of these 
methods to a 24 hour collection of the AIS signals on a selected set of ships. 

Before discussing prediction, it is necessary to describe some of the characteristics 
of the values that are to be predicted. These values are reported as discrete values. For 
some values, such as position, the underlying value is continuous, but the data itself is a 
quantized version of that continuous value. For other fields, such as the AIS field msg_id, 
the underlying value as well as the measurement are discrete. Different predictors are used 
for each of these underlying models. 

For each field, the process noise describes the deviation of the value predicted by 
the model from the actual value. The measurement noise describes the error in measuring 
the actual value. The auto- and cross-correlation functions of the model (process) noise or 
the measurement noise are not known; they are estimated. 

A third and very important point is that the error term is not the difference between 
the true state and the estimate. The error term is the difference between the estimate and 
the reported, or quantized state. This is because the estimate is used to help decide what 
was transmitted, and it is the quantized values that are transmitted. It is also interesting 
to note that if the underlying noise in the system is significantly less than the quantization 
error, then the error in the estimate can not accurately be modeled as a stationary random 
process. Depending on how the true distribution aligns with the quantization levels, the 
PMEF of the quantized value could be quite different from that of the true distribution. To 
see this, consider a quantized random walk. When the underlying true value is close to 
a quantization level, the output will likely be constant. When the true value is exactly 
between quantization levels, the output will alternate between the two levels. 

Prediction is a tool used in this work, but not the focus of the work. The devel- 


opment below is intended to introduce the concepts to support their later use as tools. 


AT 


The objective of this chapter is not to develop the “best” predictor. The objective is to 
demonstrate that the prediction assumed in the rest of the dissertation is possible, and to 


quantify the performance of the predictor. 


A. FIRST ORDER PREDICTOR 


The least complicated of the predictors considered is the first order predictor. The 

first order predictor estimates the value at time n+/ as a factor times the value at time n 
[90]: 

fin +] =aof[nl (4.1) 


For many applications, the first order predictor with ag = 1 may be the most appropriate 
predictor. Consider a remote outdoor temperature sensor with a resolution of 1/2 degree 
Celsius reporting once per second. A practical predictor may suggest future values are 
identical to current values. The first order predictor is appropriate when the process is 
modeled as discrete Markov, and the probability of transition out of the current state is 
less than the probability of remaining in the current state. This occurs when the diagonal 
elements of the transition matrix are greater than any other element in a particular row, 
Ai; > Ai; Vj #1. The AIS Navigation Status field is an example of a process that 
can be accurately modeled as a discrete Markov process with this property. The primary 
advantage of the first order predictor, and the reason it is used here, is simplicity. 

Most of the fields in the type 1, 2, or 3 AIS messages are modeled as discrete 
Markov chains with the property described above. For the type 1, 2 or 3 AIS packets a 
first order predictor is used for the following fields: Repeat Indicator, User ID, Navigation 


Status, Rate of Turn, Speed Over Ground, Position Accuracy, Course over Ground, and 
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True Heading. The residual entropy of the error from a first order predictor for these AIS 
type 1, 2 and 3 fields are given in Figure 3.8 for an example 24 hour data set. 


B. SECOND ORDER PREDICTOR 


The second order predictor is another very simple predictor. It predicts the next 


value based on the current value and a previous value using linear interpolation [90]. 





fin+a] = fin) + atin 9 (4.2) 


The second order predictor is also very simple to implement. The real problem with the 
second order predictor is that noise on the measurement can be magnified by the //k factor. 
Because of this limitation, the second order predictor is not used for any of the AIS fields. 

A better predictor is one that somehow weights new measurements based on some 
knowledge of how much noise is in the system. Predictors based on the Kalman filter have 


this characteristic. The Kalman filter is introduced in the following section. 


C. KALMAN FILTER 


The Kalman Filter is a recursive filter that provides an optimum, MMSE estimate 
of an unobserved state, provided that the system can be modeled as a Markov chain with 
linear operators relating the next state to the current state, perturbed by white noise [91]. 
In this work, the Kalman filter is used as a predictor to estimate AIS fields that nominally 
change at a constant rate, such as latitude and longitude. The Kalman filter is described 
below using an example from [91] and [92]. 

Consider predicting the one-dimensional position of a train on a straight track. The 
unobserved state consists of position and velocity. The system is modeled as ideal; having 


no external forces leads to a state transition equation of 
Lpy1 = Ope, + We, (4.3) 


where ®; is a 2 x 2 state transition matrix, k is the time index, and a is the column state 


vector with first and second elements of position and velocity, respectively. The process 
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noise vector, W,, is a zero-mean multivariate Gaussian random variable with covariance 
matrix Q, = E{w;,w}}. The vector w, is uncorrelated in time, E{w,w? } = 0, Vi £ k. 
The noise vector w, represents the deviation from the ideal model due to unmodeled forces. 
For the train example the deviation from the model could be caused by a slight incline or 
uneveness of the track, the conductor accelerating, etc. All of these effects are modeled 
with the noise term wy. 


For this example, the transition matrix is: 


1 At 
0 1 


— 


(4.4) 








Where At is the time step between successive measurements. Position at time k +1 is the 
position at time k plus At times the velocity term, and the velocity at time k + 1 is equal 
to velocity at time k. 

The observation at time k is denoted by the vector z,. The mapping of the true 


value of the state x, to our measurement Z, is through the matrix H 
Zp = Aya, + vy, (4.5) 


where the measurement noise vector v;, is a zero-mean multivariate Gaussian random 
variable with covariance matrix Ry = E{v;,v/}. As in the case of the process noise the 
vector v; is uncorrelated in time, E{v,v)} =0, Vi Fk. 

For this example assume position is the only measurable quantity. The measurement 


vector z, is then a scalar, and the matrix H reduces to a row vector. 
H= E 0| (4.6) 


The algorithm starts with an initial state estimate, %, and covariance matrix, Po, 
of the initial estimate, where P, = €{(#, — x,)(%, —2,)'}. The algorithm then follows a 
cycle of predicting ahead using the transition matrix followed by updating the prediction 
with a measurement. 

The estimate of the next state, prior to measurement, is found using the previous 


estimate of the state and the transition matrix to project ahead: 


Lesilk = Pr€uje, (4.7) 
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where the notation £,), indicates the estimate & at time j based on all the information up 
to time k. 


The estimate of the error covariance is also projected ahead using: 
Prarie = ®pP yp, P, + Q, (4.8) 
The estimate is then updated with the current measurement as: 
elk = Beer + Ke (Ze — We ®uye-1) (4.9) 


Where K,, is the Kalman gain (described below) and the term in the parentheses, z, — 
AL;,2j\n-1, represents the new information present in the measurement. 


The estimate of the error covariance is also updated with the new gain 
Py, = (LF - Be Ax) Pri (4.10) 


where I is the identity matrix. 

The Kalman gain, Kx, is selected to minimize the mean square error €{(a, — 
Brjk—1)! (Le — &p\n-1)}. The mean square error for the estimate is the trace of the matrix 
P,,. It can be shown [91], that the Kalman gain that minimizes the mean square error is 
given by: 

Ky = Pyx-1 Ay, (AP yj, + Re) (4.11) 


Inspection of Equation (4.11) reveals that when the elements of R, are large, indi- 
cating large uncertainty in the measurements, the elements of the Kalman gain are small in 
magnitude. This effectively lowers the weight of new measurements in updating the state 
in Equation (4.9) and updating the estimate of the error covariance in Equation (4.10). 
The next state would then be determined almost entirely by the previous state and the 
state transition matrix as given in Equation (4.7). A trade-off exists in selecting Q, and 
R,, between convergence and how well the filter attenuates noise. If the filter is modeled 
with too little process noise Q, or too large measurement noise Ry, the filter will take a 
very long time to converge, or may never converge, but will do a very good job of attenu- 
ating noise in the measurements. On the other hand, if the filter is modeled with too large 
process noise Q,, or too little measurement noise Rx, the filter will converge rapidly but 


do a poor job of attenuating measurement noise. 


51 


Table 4.1: Prediction algorithm based on the Kalman filter 
Xo = initial estimate 
Po» < initial estimate 


loop 
K,<= Py\,-1 Hj, (AP yp Aj, + Ti {Compute Gain} 
Pie = Lana + Ke(Ze — Ai®ne—i) {Update State Estimate} 
Pi. = UT -— KA x) Pei {Update Estimate Error Covariance} 
Lrsifk = Pp€rjx {Project ahead to next state} 
Prk = ©, Py. P;, +Q, {Project ahead Estimate of error covariance} 
end loop 


If the parameters for the Kalman filter are known, the Kalman filter provides the 
optimum MMSE estimate of the state [91]. The steps in the Kalman filter algorithm are 


summarized in Table 4.1. 


D. SPECIAL NETWORK FLAG PREDICTOR 


The SOTDMA Communication State field is structured as shown in Table 4.2 and 
Table 4.3 [1]. Clearly, most of the parameters in the table are predictable once the first 
frame has been received. Some of the parameters are effectively deterministic, such as 
time. 

The SOTDMA Communication State field of an AIS message is almost entirely pre- 
dictable. Figure 4.1 is presented to give an example of typical SOTDMA Communication 
State field values as a function of time. The figures were created using data collected from 
a single ship over a 24 hour period. Of the possible values contained in the Communi- 
cation State field, only the number of received stations and the time slot offset are not 
deterministic. 

A difficult problem in determining the value of the communication state is determin- 
ing the “timeout” value. When the timeout reaches zero, a new value is selected randomly 
by the SOTDMA algorithm. The probability of correctly guessing this new value is 1/7. 
Unfortunately, this timeout value also determines the meaning of the submessage. 

The algorithm developed for predicting the Communication State field is summa- 
rized in Table 4.4. The algorithm begins with the time of the message reception in seconds 


of day (SOD), the received data (rx), and the previous values: timeout, SOD, number of 
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Figure 4.1: SOTDMA fields 
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Table 4.2: SOTDMA communication state after [1] 
Parameter Bit Position Description Parameter Modeled 
As 
Sync State 19:18 0 UTC Direct, 1 UTC indirect, random constant 
2 Station synchronized to another 
base station 
slot time-out 17:15 Specifies the number of frames re- Random constant and 
maining until a new slot is se- deterministic 
lected, 0 indicates this is the last 
transmission in this slot 
Sub Message 14:1 Submessage meaning depends on See Table 4.3 


value of the timeout as described 


in Table 4.3 





Table 4.3: SOTDMA communication state submessage after [1] 











Slot time- Submessage Description Parameter modeled as 
out 
3,0, 7 Received Indicates the number of other sta- slowly varying random 
Stations tions currently received (between walk 
0 and 16383) 
2,4,6 slot num- Slot number used for this trans- deterministic 
ber mission (between 0 and 2249) 
1 UTC Hours is indicated in bits 13:9 approximately deter- 
hour and minute is coded in bits 8:2 ministic 
minute 
0 slot offset | Indicates the relative jump to the Uniform distribution 


time slot for transmission in the 
next frame. 


around current slot 
number 





received stations, and slot number. It first generates estimates based on prior received 


values, then demodulates the received data based on these estimates. Finally, it updates 


the stored values to be used for future predictions. 


EK. 


In this section some experimental results are presented for the prediction algorithms 
developed for the AIS field data. A 24 hour data set provided by the USCG and extracted 


RESULTS FOR AIS 
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Table 4.4: Predict communication state 


«<= sod mod 60 
{predictions} 
estTimeout = (pTimeout|i] — |(sod — pSod|i])/60 + 0.5|) mod 7 
estSyncState = pS yncState 
if estTimeout = 3|5|7 then 
estSubMessage < pSlot|i]; 
else if estTimeout = 2|4|6 then 
estSubMessage = numStations; 
else if estTimeout = 1 then 
estSubMessage = |(sod mod 24*60*60)/(60*60) | *2°+|(sod mod 60*60) /60 «2? | 
else if estTimeout = 0 then 
estSubMessage <= 2248 
end if 
{demodulate} 
{timeout, syncState, sub} <= demod(r, est) 
{update for future predictions} 
pSod|i] <= sod 
ptimeout|i] = timeout 
psyncState|i] = syncState 
if timeout = 3|5|7 then 
pSlot|i| = subMessage 
else if estTimeout = 2|4|6 then 
numsStations <= subMessage 
else if estTimeout = 0 then 
newSlot < (slot{i] + subMessage — 2250) mod 2250 
plimeout <= 7 
end if 


with the AlSUser application is used to generate the results presented in this section. The 
receptions are sorted by transmitting ship using the MMSI field. The data is then passed 
to a predictor that attempts to predict the next message field-by-field given the previous 
messages. 

The receive station that generated this data received reports from 192 unique ships 
on the collection date. Results from three of the 192 ships are presented below. 

Three ships were selected for analysis. Two ships under way with two different 


reporting precisions for the latitude and longitude fields. Ship A has a reporting precision 
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100 counts and Ship C has a reporting precision of 10 counts. An additional ship, Ship B, 
is selected that is not under way; this ship has a reporting precision of 1 count. 

The length of time from a measurement to an estimate is variable. For example, if a 
ship reports every six seconds, the error will nominally be based on a six second prediction. 
In practice, the data set provided is missing some reports. This is particularly true when a 
ship is far away from the receiving station. An additional criterion in ship selection for this 
example is contiguous measurements, or a small number of missing measurements. The 
time between measurements for the first 2000 measurements of ships A, B and C is shown 
in Figure 4.2. Although there are instances of up to 90 seconds without a report, the plots 
are scaled such that time differences of greater than 40 seconds are not shown to illustrate 
the detail for the nominal differences. Notice the dip in time difference for Ship A around 
measurements 400 and 750. This represents a time when the ship is changing course, 
and the reporting rate increased accordingly. Ship A and Ship C are moving, and have 
a nominal reporting interval of six seconds. Ship B is nominally stationary, and has a 
nominal reporting interval of ten seconds. The following conclusions can be drawn from 
Figure 4.2: there is significant variation around the nominal reporting interval, and the 
data set is likely missing some reports. 

Let the error, e{k], be defined as the difference between the reported (quantized) 
value, q(a,) at time k and the estimation of the reported value, q(#j\x-1), at time k based 


on all measurements before time k: 
elk] = q(@xjn-1) — (Lr) (4.12) 


The q(-) operator represents the quantization operation resulting in the reported quantized 
value rather than the actual unquantized value. 
The entropy (in bits) of the error term, e, with discrete probability distribution 
(pi, P2,---,PK) is defined as 
K 
H(p) = So pk log,(1/px) (4.13) 


k=1 
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Figure 4.2: Time Between Position Reports 


where we define 0 log,(1/0) = 0. The distribution of the error term can be estimated with 


a normalized histogram. 


1. Longitude Field 
The latitude and longitude fields are unique in that there are many bits used to 


represent the value, yet there is little information contained in additional measurements 
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after one measurement has been received. Prediction results for the longitude field are 
presented below. The results for the latitude field are similar. 

A two-state Kalman filter is used to predict the longitude field. The state transition 
matrix is identical to the example of predicting the location of a train on a straight track 


discussed in Section IV.C 
1 At 


Os 


k= 


(4.14) 








where the time step At is taken to be one second. 

The process noise must be modeled large enough to model maneuvers that a ship 
may perform; this is effectively modeled as noise in the velocity term. Consider an un- 
known acceleration term, a, modeled as a Gaussian random variable. The w, term in 


Equation (4.3) becomes 
Ae 
aK (4.15) 








The random acceleration, az, is modeled as a zero-mean Gaussian random variable with 
standard deviation o,,. For this experiment o,, is chosen to be 0.71 counts/sec?, where 
the units are binary counts, i.e. the raw values that are transmitted. There is no need to 
convert to physical units. 


The covariance matrix for the process noise Q, is then 








L/4 Age. 1/9AP 
Q,, = 0.5 / / (4.16) 
L2ae Ae 
The variance of the measurement error is set to 
Ry = 20; (4.17) 


with units of counts squared. This value is determined empirically by looking at many 
sets of data. The variance of the measurement error represents both errors from the ship 
estimate of its position and quantization errors. 

The initial estimate Zp has the first element set to the initial measurement with the 


second element, representing velocity, set to zero. The initial error covariance is set to 


P, - ie | | (4.18) 
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representing a large uncertainty for the initial estimate of the velocity term, and a relatively 
small uncertainty for the initial estimate of the position term. 

The prediction results for the longitude field are presented for each of the three 
ships in Figures 4.3, 4.4, and 4.5. For each of the ships, the estimates of longitude and 
longitude rate are plotted. The error in the longitude estimate is also plotted versus time. 
A histogram of this error with a measure of the entropy of the error is also presented. 

Figure 4.3 shows the prediction performance of the Kalman filter based predictor 
for a ship traveling at a nominally constant rate with a reporting precision of 100 counts. 
Figure 4.3b illustrates two clear changes in longitude rate. Notice in Figure 4.3c and 
Figure 4.3d that the distribution of values around the estimate is effectively just three 
values, the nominal estimate and one on each side. Also notice that the rate changes that 
occurs at approximately 2200 and 4600 seconds from the first measurement do not cause 
an increase in error in the estimate. The ability to track rate changes is incorporated 
into the process noise matrix @,. The entropy of the error in the estimate, or the added 


information from an additional measurement is approximately 1.2 bits. 
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Figure 4.3: Longitude Field Prediction Performance Ship A 
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Figure 4.4 shows the results from a ship that is not moving that reported with a 
precision of 1 count. As indicated in Figure 4.4b, the velocity of Ship B is effectively zero 
for this data set. Notice in Figure 4.4d that the distribution of values around the estimate 
is effectively five values; two on each side of the estimate. The entropy of the error of the 
estimate, or the added information from an additional measurement, is approximately 2.1 


bits. Although the error in the estimate is lower for Ship B than for Ship A, because of 


the higher reporting precision, the entropy of the error is actually higher. 
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Figure 4.4: Longitude Field Prediction Performance Ship B 


Figure 4.5 shows the results of a ship with similar motion to that of Ship A but 
with a reporting precision of 10 counts. In this example the error in the estimate increases 
as the magnitude of the rate of change in the longitude value increases. Changes in the 
velocity term of the Kalman filter are modeled with a noise term. It will take some time 
for the Kalman filter to settle on a new estimate for velocity; during this time the error 
in the estimate will be of a larger magnitude. There is some unexplained non-uniformity 
in the distribution of the error seen as some banding, or a set of outliers, in the error plot 

that has not been accounted for. This non-uniformity can also be seen in the histogram 
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at around +50 counts and —50 counts as two small humps. Notice in Figure 4.5d that the 
distribution of values around the estimate is effectively 13 values, six on each side of the 
estimate. For this case the entropy of the error estimate, or the added information from 


an additional measurement is approximately 3.5 bits. 
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Figure 4.5: Longitude Field Prediction Performance Ship C 


Some observations follow: Both the reporting precision and whether or not a ship 
is moving have an effect on the entropy of the AIS fields. It is apparent from the data 
from Ship C that the error can not necessarily be modeled as a stationary random process. 
The distributions for the longitude field will theoretically contain 27° elements, but for all 
practical purposes only a small number of elements are required to represent the distribu- 


tion when prediction strategies such as those presented here are used. The significance of 
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this will become more apparent in Chapter VI which discusses the complexity of a MAP 


detector. 


2. SOTDMA Communication State Field 

The SOTDMA field is predicted using the algorithm summarized in Table 4.4. The 
residual entropy in the error from prediction of the Communication State field is not 
dependent on the reporting precision or the reporting rate. Therefore, only the results for 
one of the ships (Ship B) are presented. 

The results from the Sync State field are trivial. This field did not change for any 
of the three ships for the duration of the data sets. The calculated entropy for this field 
would therefore be zero. 

The prediction error for the Timeout field is presented in Figure 4.6a with a his- 
togram showing the distribution of the error in Figure 4.6b. The Timeout field was cor- 
rectly predicted approximately 79% of the time. For this data set the timeout was zero 
approximately 18% of the time. When the timeout is zero, the next value of the timeout 
is randomly selected from the set {1,2,3,...,7}. Therefore the probability of an error in 
prediction is 6/7. There are two conditions under which timeout is incorrectly estimated: 
the slot timed out and a new timeout is randomly selected, or a message is missed such 
that the previous timeout is unknown. When the slot times out, a timeout value of seven 
is arbitrarily selected. For this reason, most of the errors in the estimate are on one side 
of the distribution. 

The results for the communication state submessage is shown in Figure 4.7. Close 
examination of Figure 4.7a reveals three types of errors corresponding to errors estimat- 
ing the time of day, the slot, or the number of received stations. If the timeout value 
is incorrectly determined, that incorrect estimate will be used for determining which of 
the three potential values of the submessage to estimate, leading to three distinct errors. 
Observe from Figure 4.7b and 4.7c that the distribution of values around the estimate is 
tight for approximately 75% of the values with the remaining 25% of the values are spread 
over a large number of count values. There is a somewhat uniform distribution of errors 
among the 2250 possible slot values, another Gaussian type distribution for the number 
of receiving stations and then a time-varying distribution around the hour and minutes 
of the Time field. The calculated entropy for the submessage is approximately 4.4 bits, 
significantly less than the 14 bits that are transmitted for this field. 
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Figure 4.7: Communications State Submessage Field Prediction Performance Ship B 
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Figure 4.6: Communications State Timeout 
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F. CHAPTER SUMMARY 


This chapter develops an efficient method to predict the AIS data fields using three 
different types of predictors. An algorithm is developed for prediction of the AIS Communi- 
cation State field. The results of these predictors are then presented on three representative 
data sets. The next Chapter discusses methods of incorporating knowledge of bit-stuffing 
into a MLSD. In the following Chapter Maximum a posteriori detection is discussed; here, 
the importance of distribution knowledge on performance and the importance of distribu- 


tion reduction on computation time is illustrated. 
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V. MLSD FOR BIT-STUFFING 


This chapter investigates maximum likelihood sequence detection (MLSD) for a bit- 
stuffed data source. The AIS uses the HDLC protocol which has bit-stuffing as one of its 
components. 

Bit-stuffing is used in various communications protocols, typically as a way to re- 
move ambiguity in the detection of control flags or as a method to aid clock recovery. Typ- 
ically a demodulator makes hard symbol-by-symbol decisions without incorporating the 
knowledge that the data source is bit-stuffed. A demodulator making symbol-by-symbol 
decisions is clearly not optimum for a source that has been bit-stuffed; symbol-by-symbol 
decisions ignore the knowledge that the source has been encoded. In this chapter an opti- 
mum approach is investigated. The maximum likelihood sequence detector for a bit-stuffed 
data source is developed; starting by modeling the transmitted signal as a Markov chain. 
Approximate performance bounds for the MLSD are calculated. Finally, a trellis receiver 
structure is developed, and perform simulations using the Viterbi decoder algorithm are 


performed. 


A. INTRODUCTION AND MODEL 


Bit-stuffing is the procedure of introducing an extra bit after a fixed number of 
consecutive ones or zeros [93]. This chapter is concerned with the insertion of a zero 
after k consecutive ones; let this be indicated as 1* bit-stuffing. The High-level Data Link 
Control (HDLC) protocol is an application where 1° bit-stuffing is used. Bit-stuffing has 
two obvious benefits: it aids in clock recovery, and it removes any ambiguity in locating 
the start and stop flags, 0x7E (01°0), for HDLC because these flags are not bit-stuffed. 

Typically unstuffing is handled after demodulation. The demodulator is designed 
as if the data bits were independent, making symbol-by-symbol decisions. A demodulator 
making symbol-by-symbol decisions is clearly not optimum for a bit-stuffed source. 

An alternate approach in demodulation is to take advantage of the knowledge avail- 
able for the signal. For example, in 1° bit-stuffing the sequence 111111 is not a valid 
sequence; if a demodulator outputs that sequence it would be known that an error oc- 
curred. A maximum likelihood sequence detector will appropriately take advantage of this 


information; this detector is efficiently implemented using the Viterbi algorithm. 
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It should be noted that once the demodulator has made hard decisions nothing can 
be done to improve the error performance. For example, with 1° stuffing if six ones are 
received in a row, there has been an error, but the error could be in any of the six bits. 
There is a 1/6 probability of correcting the error, but there is a 5/6 probability of creating 
another error. Another approach is to use soft outputs to help decide which of the six 
possible bits were in error; this technique is not optimal, and will not be discussed here. 

Bit-stuffing is similar to run-length-limited (RLL) codes. A (d,k) run-length-limited 
code has the property that consecutive 1’s are separated by at least d but no more than 
k; zeros. Bit-stuffed data has the property that consecutive zeros are separated by at least 
zero ones, and no more than & ones. Many authors have investigated error control coding 
for bit-stuffed data or RLL codes [94,95]. We [96] are not aware of any other attempts to 


examine the performance of a maximum likelihood sequence detector for bit-stuffed data. 


1. Bit-Stuffing 

Let a be a vector of information bits. A new stuffed vector a, is created by bit- 
stuffing the original vector a. For this dissertation, bit-stuffing consists of the insertion of 
a zero after each consecutive string of k ones, where the leftmost bit is transmitted first. 


Let B;, denote the bit-stuff operator. For example, if 
a = 011101, (5.1) 


then 
a, = Bo(a) = 0110101, (5.2) 


where the overbar indicates an inserted bit. Let B,' denote the bit unstuffing operator 
such that 
a= B,"(a,). (5.3) 


The received vector of decision statistics can be defined as 


r=a,+w, (5.4) 
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Figure 5.1: Model for 1° Markov Process Data Source 


where w is a vector of samples from a bandlimited white noise process. The samples in w 


have variance o?. 


2. Modeling the Bit-Stuffing Process 

The sequence a, can be modeled as a Markov chain. A Markov model for a, = 
B;(a), where a contains independent equally likely symbols, drawn from {0,1}, is shown in 
Figure 5.1. Note that for each state with the exception of state S; the probability of a one 
is equal to the probability of a zero, Pr(0) = Pr(1) = 0.5. For state Ss, a zero is transmitted 
with probability one. It is this extra piece of information that a demodulator can exploit 


to obtain a slightly lower bit error rate (BER) than that of the symbol-by-symbol detector. 


A Demodulators 
The first demodulator examined is a symbol-by-symbol maximum likelihood detec- 


tor. This detector selects a,[k] based on the following criteria 


a.[k] = arg max_Pr(r|k]la,[A]). (5.5) 


as€{0,1} 


It should be clear that this first detector is not able to incorporate knowledge of valid or 


invalid sequences; decisions are based only on the current sample. 


67 





So Os 





Sy 2 
S» ° 
Ss ° 
Sy ° 
S5 ° 


Figure 5.2: Trellis Structure for 1° Bit-Stuffed Markov Process 


The second demodulator considered is the maximum likelihood sequence detector. 


This detector selects the entire vector a, based on the following 


a, = arg max Pr(rja;) (5.6) 


sCAs 
where A, is the set of valid stuffed sequences. The MLSD incorporates what is known 
about the data source into the demodulator. The MLSD is efficiently implemented using 
the Viterbi algorithm. An example of a trellis for decoding 1° bit-stuffed data is shown in 
Figure 5.2, where the process is in state S,, if the previous n bits were ones. Transitions 
due to a “1” are indicated with a dashed line, and transitions due to a “0” are indicated 
with a solid line. The only state with more than one branch entering is state So. The add- 
compare-select function in the Viterbi algorithm needs only to be applied at this state, 


reducing the computational complexity of the algorithm. 


B. CATASTROPHIC ERRORS 


Bit-stuffing is a catastrophic encoding process; a finite number of errors in the 
encoded stream a, can lead to an infinite number of errors in the recovered data stream 
a. Most applications that use bit-stuffing are packet-based so an infinite number of errors 


is reduced to a single packet error. An interesting question is “What is the probability of 
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a catastrophic error given a single error in the encoded stream a,?” A catastrophic error 
occurs in the unstuffing operation when a bit is either removed erroneously or not removed 
when it should be. 

There are two first order phenomena where a single flipped bit leads to a catastrophic 
error. First consider a bit removed erroneously due to a single error. This requires a special 


sequence. For example, consider 1° bit stuffing and the sequence below: 
011011010110011110011111010 


Any error in the bits with the overbar will cause a bit to be erroneously removed leading to 
a catastrophic error in the unstuffed output. Let these bits with the overbar be called type 1 
bits. A type 1 bit for 1* stuffing is a zero that can be flipped to create a non-overlapping 
sequence of k sequential ones. 


Next, consider a stuffed bit erroneously not removed from the same sequence. 
011011010110011110011111010. 


Errors in the bits with the overbar will cause a stuffed bit to be erroneously left in the 
stream leading to a catastrophic error in the unstuffed output. Let these be called type 2 
bits. A type 2 bit for 1’ stuffing is any of the k ones in a sequence of k sequential ones. 
The probability (P.,-) of a catastrophic error beginning at a particular bit is ap- 
proximated by 
Hee fpaliec Papel; (5.7) 


where P,; and Py are the probabilities that a particular bit is a type 1 or type 2 bit, 
and Py. is the probability of any particular stuffed bit is received in error. This is an 
approximation and not an equality because higher order error events have been ignored, 
such as two errors leading to the erroneous removal of a bit, or the case of a second error 
canceling the effects of the first. 


The probability of a frame containing a catastrophic error is 
Pee = (Pe) (5.8) 


where n is the length of the frame. 
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Figure 5.3: Pe and P.fe versus Frame Length n, 1° Stuffing and P.». = 1074 


A plot of the probability that a catastrophic error will begin at a certain bit, P.e, 
along with the probability of a catastrophic error occurring during a frame, P.fe, are shown 
in Figure 5.3. Notice that P.. grows rapidly with frame length and then quickly levels out 
to be approximately constant. For small probability of stuffed bit error and frame lengths 
greater than around 20, P.fe will increase linearly until the probability of multiple errors 


becomes significant at which point it will slowly flatten out and approach one. 


C. CALCULATION OF Py 


There has been some interest in computing P:2 by those interested in calculating the 
probability that there will be N or more stuffed bits in a given sequence. Nolte, et al., [97] 
performed brute force calculations of the number of non-overlapping runs of k ones in n 
bits for their work in latencies in the Controller Area Network (CAN) bus. An alternative 
approach is demonstrated below. 

Let Np. be the number of non-overlapping runs of & ones in n bits. Begin by 
calculating the probability of one or more non-overlapping runs of k ones in a stream of n 
random bits Pr(N,,, > 1). 

Markov chains can be used for this calculation without resorting to brute force or 


Monte Carlo approaches [98]. Consider a sequence of n independent bits where Pr(1) = p, 
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Figure 5.4: Markov Chain for Calculating Pr(V,, > 1) 


and Pr(0) = 1 — p. A state diagram for the Markov chain used for finding Pr(N,4. > 1) 
is shown in Figure 5.4. The states are defined such that the Markov chain is in state S;, 
when N,,, > 1 and in state 5;, for 0 <7 < k, when N,,;, = 0 and the last 7 bits are ones. 
Except for state S;, each of the states advances with probability p or returns to state So 
with probability 1 — p. State S; is referred to as the absorbing state; once entered it is 
never exited. 

The (k + 1) x (k +1) state transition matrix for the Markov chain is 


Le pO pe see dee 0 
p 0 0 aie 0 O 
0 p 0 0 O 
Il, = ; ; ; ’ (5.9) 
0 0 
0 0 p 0 O 
0 0 0 0 p. ~s 


where II; {2, 7] = Pr(Si41 = 2|S; = 7). The initial state probability vector is €) = eo, where 
e; is a length (k + 1) unit column vector with a 1 in the (i + 1)-th position and zeros 
elsewhere, €9 = E 00... 0 0. - The state probability vector after n bits is given 
by 

E,, = W1"eo, (5.10) 
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Figure 5.5: Markov Chain for Calculating Pr(V,4 > L) 


where €,, is the state probability vector at time n. The probability of at least one sequence 


of k ones in n bits is the probability of being in state k, which can be written as 
Pr(Nnp > 1) = ef Teo. (5.11) 
The probability that there are no sequences of & ones in n bits is given by 
REG = Or Prag = 1). (5.12) 


The above can be extended to determine Pr(N,,, > L), the probability of L se- 
quences of & ones in a stream of n random bits. A state diagram for the Markov chain 
used for finding Pr(N,,, > L) is shown in Figure 5.5. The extended Markov chain consists 
of L cascaded rings, L — 1 rings of identical form ending with a ring of the form identical to 
Figure 5.4. The states are defined such that the Markov chain is in state S;_1, if Nny, > L, 
in state S;;, for0O <7 < L—1if N,, =i+1 and the most recent k bits are ones, and in 
state S;; forO0 <i<L—landO<j <kif Ni, =7 and the last j bits are ones. 

The L(k +1) x L(&k +1) transition matrix for the extended Markov chain is 


Il,’ Z Z 
N Z 
W,=|Z N Z (5.13) 
Il,’ Z 
Z N 


2 


where IT,[j+(k+1)i,m+(k+1)]] = Pr{.S(¢+1) = S;;|S(t) = Sim)}, Z is a (kK+1) x (k+1) 
zero matrix and N is a (k+ 1) x (k +1) zero matrix with the exception of two terms in 


the upper right corner representing the transition probabilities to the 5,9 and Sj; states. 


0 0: cL 
0 Os - 69 
N= (5.14) 
0 0 
0 0 0 O 
II,’ is identical to the Ty matrix without the absorbing state Sin: 
l=p. b= po l= ee bap. 0 
p 0 > Sher SO 0) 
0 0 a vOo- VO 
II,’ = i (5.15) 
' 0 O 
0 0 0 p 0 0 
0 0 0 0 pO 
The probability of being in a particular state after n bits is given by 
E,, = It" eo, (5.16) 
and the probability of at least DL sequences of & ones in n bits is then 
Pr(Nng 2 L) = €7(441)-1114"e0, (5.17) 


as before. 

Consider determining P)2, Piz is the probability that a specific bit in a, is one of the 
particular bits that when flipped will cause a “type 2” catastrophic error. This is found 
by summing over the probability of exactly | sequences of length k, multiplying by the 
number of bits (k) in those sequences and dividing by the total number of bits (n). Let B 
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describe the event that a bit is in arun k ones. Py is then 


[n/k| 
Po = >> Pr(BlNne =!) Pr(Nne =!) 
l=1 
yp n/t 
es S> Pr(Nnk = DL. (5.18) 
l=1 


Pr(N,,,~ =!) can be calculated using (5.17) as follows 
Pr(Nn~ = 1) = Pr(Nne > 1) — Pr(Nn, > 1+ 1). (5.19) 


It should be noted that Philippou and Makri [99], and independently in the same 
year Hirano [100], developed a formula for calculating Pr(N,~. = 1). Unfortunately it is 


not clear how to extend their results to find Py. 


D. CALCULATION OF Pi, 


Consider the calculation of P;,;, P;, is the probability that any bit in a, is one of the 
particular zero bits that when flipped will create a non-overlapping sequence of k sequential 
ones leading to a “type 1” catastrophic error. For example, for 1? stuffing the sequence 
0100011 has type 1 bits located at the zeros with the overbar, but the zero with the hat 
is not a type 1 bit, from the non-overlapping requirement. Calculating the probabilities of 
these locations is unfortunately more difficult than the case of type 2 errors, but the same 
general technique that was used to determine Pz can be applied. The sequences take the 
form 1/-!01*-4, where j is the location of the zero in the sequence of ones, 7 = 1 being the 
leftmost bit position. Let Ne ‘Ol be the number of runs of 17-101" that do not overlap 
with a run of k ones. A Markov chain for finding Prine, oP Se Te I eh) 
for 1* bit-stuffing is shown in Figure 5.6. The Markov chain can be extended as before to 


calculate the individual terms, and find 


2x [n/k] . . 
Ris PN. SOE (5.20) 
l=1 


1 
n 
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Figure 5.6: Markov Chain for Calculating ee Se ED cob) tor de 
Bit-Stuffing 


Finally, Equation (5.7) may be used to find P».. To understand the limits on the sum in 
Equation 5.20 consider the following sequence for B?: “010101010” with n = 7 and k = 2. 
The number of type 1 bits in this sequence is four. The limit represents the maximum 


number of special bits that could occur in a sequence. 


E. BOUNDS ON BIT ERROR RATE PERFORMANCE 


Deriving an analytical expression for error performance is slightly more complicated 
than finding an error performance bound for a linear convolutional code. Assuming the 
all zeros sequence is sent results in an overly pessimistic error bound. This is because the 


invalid codes are relatively far in Euclidean distance from the all zero sequence. Forney [51] 
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showed that a tight lower bound for probability of error for a MLSD is 


Fhe ee KoQ(dmin/2e), (5.21) 


where o is the standard deviation of the noise component of the decision statistic. For 
this example, Ko is the probability that an error of Euclidean distance d,,;n from the given 
input sequence will lead to a valid sequence. A lower bound on performance can be found 
relatively easily using the results from Section V.B. 

The encoding format from bit-stuffing reduces the set of possible valid sequences. 
To calculate a lower bound, first find the sequences that have a distance of dmjn to one 
of the invalid sequences. The lower bound on probability of bit error is then the average 
number of din error events divided by the number of bits. These events are similar to 
the “type 1” events, but not identical. For example, for 1° bit-stuffing the sequences that 
have distance of 1 from an invalid sequence are 111110, 111101, 111011, 110111, 101111 
or 011111. The events of interest are the zeros in the 1° stream a,, where the number of 
consecutive ones to the left and the right sum to at least k = 5. Let Z,,;, be the number of 
zeros in a sequence of n bits stuffed using B; with k ones to the left and right of each zero. 
As before Pr(Z,,, > 1) is found using Markov chains. This involves a Markov chain very 
similar to the one used to find “type 1” errors. Figure 5.7 illustrates this for 1° stuffing. 
Note that while this is very similar the Markov model used in calculating “type 1” errors for 
1* stuffing, it is distinct beyond the change in bit-stuffing level; there is no non-overlapping 
requirement. In addition, after k ones a zero is inserted, so the chain jumps to a different 
state. It is important to notice that the Markov chain is run on data prior to bit-stuffing 
(a), but it is predicting values for the bit-stuffed data (a,). The sequence of bits where 
this comes into play is 1* or 01*. In both of these cases it is known from the bit-stuffing 
rule that the next bit in a, will be a zero and it will be one of the special zeros. Similar 


to P the results from the Markov chain are used to calculate Ko 


2x[n/k] 
Kala SP Zeal (5.22) 


l=1 
Figure 5.8 shows calculated values of Ko for a few values of k versus the number 
of bits n. Clearly smaller values of k and larger values of n will lead to a larger potential 


performance improvement for the MLSD over symbol-by-symbol detection. 
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Figure 5.8: Ko versus n 
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Figure 5.9: Illustration of Stuffing, Demodulation and Unstuffing Process 


The lower bound (5.21) is tight for high SNR. For low SNR the neglected terms 
representing higher order events, more than one bit in error, become quite large and it is 


no longer a useful bound. 


F. SIMULATION 


1. Results 

Data generated according to (5.4) is used as input to two demodulators: the symbol- 
by-symbol maximum likelihood detector and the maximum likelihood sequence detector 
described above. Each demodulator generates an estimate, @,, of the vector a,. The vector 
a, an estimate of a, is then generated using (5.3). See Figure 5.9 for an illustration. 

Figure 5.10 shows the BER performance of the two demodulators for 1° bit-stuffed 
data with frame length n = 168 prior to unstuffing (a,) and after unstuffing (a). Fig- 
ure 5.11 shows the BER performance of the two demodulators for 1? stuffing. For 1° 
bit-stuffing and a frame length of 168 the term Ko is approximately 0.9458. For frame 
lengths > 20, the BER of a, using MLSD is approximately independent of frame length. 
Clearly from Figure 5.8 and (5.21) the lower bound on BER for a, using MLSD is a function 
of frame length for small values of n. 

For reasonably high SNR the BER performance for a, using MLSD is close to the 
lower bound Ko Q(,/2E;/No). This represents a modest gain in BER for the stuffed vector 
as. 

Figure 5.12 shows the gain in performance as a reduction in required E;/No to 
achieve an equal BER. Gains using the MLSD are on the order of 0.05 dB for a,. For a 
the gains using the MLSD are from 0.2 — 1.0 dB. 
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Figure 5.10: BER Curves Comparing MLSD and Symbol-by-Symbol Decisions for 1° Stuff- 
ing 


2. Interpretation and Discussion of Results 

The results from the previous section are BER curves for the stuffed vector a,. It is 
interesting to compare the performance of the detectors with the unstuffed vector a. The 
number of errors caused by a catastrophic error is a function of frame length. A longer 
frame leads to a larger penalty for a catastrophic error. For this simulation a frame length 
of 168 bits is used. 

The BER performance difference between the MLSD and the symbol-by-symbol 
detector is greater when measured for a then for a,. As described above, a single error 
in a, can lead to multiple errors in a. This is the reason for the higher BER for a. The 


vector a contains the original information so it is generally the BER curve of @ that is 
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Figure 5.11: BER Curves Comparing MLSD and Symbol-by-Symbol Decisions for 1? Stuff- 
ing 


important. The reason for the performance improvement is that the MLSD will decrease 
the likelihood of some of the decisions that would lead to catastrophic errors. Consider 1° 
stuffing and the sequence a, = 101111010. Symbol-by-symbol decision errors at either of 
the two over bar zero positions would lead to a catastrophic error in a. An error event of 


2dmin Would be required of the MLSD in order for a catastrophic error to take place. 


G. CHAPTER SUMMARY 


This chapter investigates MLSD for a bit-stuffed data source. The superior BER 
performance of MLSD versus symbol-by-symbol decisions is demonstrated, and a lower 
bound on performance is derived. The BER performance for both the stuffed stream a, 


and the unstuffed stream a have been investigated demonstrating the effect of catastrophic 
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Figure 5.12: Gain (dB) Comparison of MLSD and Symbol-by-Symbol Decisions for 1° 
stuffing 


errors. The probability of a catastrophic error for a given stuffed bit error rate Pye is 
derived. A bound for the BER performance of the stuffed stream a, is developed by 
calculating probability of sequences with a Hamming distance of 1 from an invalid sequence. 

For this analysis the case where the number of stuffed bits is known is not considered. 
It is possible that the number of stuffed bits for a particular frame would be known. For 
example, a protocol such as HDLC has a fixed unstuffed frame length while the stuffed 
frame length for HDLC can be calculated by counting the bits between the start and stop 
flags. The difference between these lengths will yield the number of stuffed bits. It is the 
author’s belief that this additional knowledge can be used to further increase performance, 
particularly in the unstuffed vector a. 

The next Chapter investigates the multicategory classification problem of Field- 
Based MAP detection. 
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VI. FIELD-BASED MAP DETECTION 


This chapter develops the MAP detector for an individual field of symbols. Here, a 
field is defined as a consecutive sequence of related symbols. A typical transmitted signal 
may contain many consecutive fields; here we consider the reception of a single isolated 
field through an AWGN channel with no ISI introduced at the transmitter. Later in the 
chapter the special case of reception of multiple synchronous signals is investigated. Both 
of these cases result in a multicategory classification problem for a multivariate Gaussian 
density function. 

There are applications where data are transmitted as a sequence of fields, with 
each field containing a number of symbols. Here, the field values are independent, but 
the individual symbols contained in a field are not independent. For the example of AIS, 
the a priori knowledge available is for a field of symbols rather than priors on individual 
symbols. A detector that takes advantage of the additional information contained in the 
field a priort probabilities will lead to a lower probability of field error than that of either 
the MAP or ML symbol-by-symbol detector. Given that the fields are independent with 
no ISI, the optimum detector for a field of symbols can be developed by observing only the 
current field. 

For this development the interest is in minimizing the probability of field error. 
While in general this will tend to also minimize probability of bit error, it will not always 
be the case, depending on the prior probabilities for particular field values. Field error is 
often a more natural measure of performance. If the field values contain the important 


information, one error or one hundred errors within a field result in the same cost. 


A. SINGLE-CHANNEL FIELD-BASED MAP DETEC- 
TOR 


Consider an information vector a, where there is a priori knowledge of the rela- 
tionship between a sequence of information symbols (a field). For example, we may know 
Pr(ag, @k41,---;@k+n = b) where b is one particular realization of the n-tuple of all possible 


information symbols. The information vector a is mapped to a symbol vector a. This 
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vector is received in the presence of AWGN resulting in the received vector 
r=a4z, (6.1) 


where the elements of z are zero-mean independent identically distributed (i.i.d.) Gaussian 
random variables. 
The maximum a posteriori (MAP) detector selects the field that maximizes Pr(aly) 
[3, 101] 
a = arg max|Pr(aly)]. (6.2) 


This decision rule results in the minimum probability of field error. 


1. Model 
This development of the multiple decision MAP field detector loosely follows ex- 
amples in [101] and [102]. Consider the reception of a length N sequence of symbols in 
AWGN 
y=ar+zZ, (6.3) 


where @ is a N x 1 vector representing a field. Each element of @ is a member of {—A, A}. 
The noise vector z is an N x 1 vector of zero-mean i.i.d. Gaussian random variables with 
variance o”. The vector a represents one of the 2% possible field values. Let only M < 2% 
field values have a non negligible probability of occurring. For each of the M likely field 
values, there is an associated prior probability of occurrence, Prla,] = py # 0. 

Each element of the vector z has the same variance, and is independent of each of 


the other elements. This leads to the conditional PDF of y given by 





Selig Dl ap 
plyla) = (2n02) "exp (<= SEY a)) (6.4) 
The MAP detector selects the particular field @ € {a@), Q2,...,@.} that maximizes 


Pr(aly). Using Bayes’ theorem and eliminating terms not dependent on a, it can be shown 
this is equivalent to selecting a that maximizes p(y|@) Pr(@). This leads to the following 


decision rule: 


& = arg max[Pr(@)p(y|a@)| 
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Where A is the set {@1, Q2,...,@.}. From Equation (6.4) this leads to: 


—(y— a)" (y - 2))] . 


202 





a = arg max [Pray exp ( (6.5) 


Any monotonic function of the right hand side of Equation (6.5) will lead to an identical 


result. The form of Equation (6.4) leads to the natural logarithm as a good choice, resulting 


yoey ya) . 
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in 





& = arg max in(Prfad) = ( (6.6) 


2. Discriminant Functions 

Another way to pose the same problem is in terms of discriminant functions. Dis- 
criminant function is a term borrowed from the pattern classification field. The discrim- 
inant function, denoted as g;(y), is a function of the observation vector that leads to the 


decision rule of selecting @ = ay if 


mY) > gly) WAR. (6.7) 


For the MAP detector 


Gi(y) = Prlaa|p(ylai); (6.8) 


this choice for discriminant function is not unique, Equation (6.7) remains valid if g;(y) 
is replaced with f(g;(y)) for any monotonically increasing function f(-). Because of the 
form of p(y|a;) given in Equation (6.4), a logical choice for f(-) is the natural logarithm, 


resulting in a new discriminant function 


Atos Pyedy (# a =) | (6.9) 
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Observe from Equation (6.9) that the discriminant functions consist of two terms: a 
squared Euclidean distance measure weighted by the noise power and the a priori proba- 
bility. If Pr[a,| = Pr[a;| Vk,1 then this results in the minimum distance decision rule of 
the ML detector. Note that unlike the ML decision rule, knowledge of the noise power 0? 
is required for the MAP detector. It is clear from Equation (6.9) that for large values of 
o” the decision rule is weighted heavily by the prior probabilities. 
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Equation (6.9) can be expanded to write the discriminant function as 


(6.10) 





y y+ aja; — 2y"a; 
20? 


sity) = In(Pr[ox)] ~ ( 
The term y’y can be removed because it is common to each discriminant function, leaving 
a new simplified discriminant: 


Taq. 
ey es (6.11) 





gi(y) = 07 In(Pr[e,]) — 


If each of the signals have equal energy, the a/ a; terms can be removed from Equa- 


tion (6.11), and the discriminant function can be further simplified to 
gi(y) = 07 In(Prlai]) + y" a4. (6.12) 


The decision region, R,;, consists of the values of y for which @, maximizes (6.5). 


The boundaries of these regions are the (V — 1)-dimensional hyperplanes where y satisfies 


guy) = gly). (6.13) 


To show that the boundaries are hyperplanes, consider the discriminant functions for the 
general MAP case. The boundary is then given by the values of y satisfying the following 


equation 


n(Pr{ox,)] — (“ = an)" (y - 0) Sie ene (“ = au) "(y — 2) 6) 
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Expand the above and cancel common terms to obtain 


afar, — 2y? a 











In(Pr[ax]) — 552 = In(Pri|a|) — 552 ‘ (6.15) 
and finally collect the y terms to obtain 
1 1 Pria 
y' (ay — a) = 5k. — sou! OH —o7 ln (Fes) ; (6.16) 
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The right hand side of Equation (6.16) is a constant scalar. The decision boundary that 
results is a hyperplane parallel to the vector y which satisfies y?(az — ay) = 0. 

If the field values are all equal energy, the discriminant functions in (6.12) can be 
used, resulting in boundaries defined by 


y" (Qi — a1) = 07 ln (sie) (6.17) 


Pria,| 

There are two special cases which guarantee reachability of all decisions for some 
specific value of y (a value of y exists such that decision @; is made Vl). The first case 
is equal energy symbols; the second is equal probability sequences. In the general case of 
unequal energy symbols and unequal probabilities, there may be some decisions that will 
never be made no matter what the value of the observation vector y. 

For the equal energy case: 

Hy, 
y! (a, — a) ZC (6.18) 
ai 
Let y = Bax, a scaler B can always be found such that decision & is made. 

For the equal probability case, the decision rule results in the minimum distance 
decision. Setting the observation vector to a, will result in a distance of zero and decision 
k; will be made. 

Figure 6.1 illustrates the decision regions for the two symbol case in R?. Notice the 
direction the decision boundary moves as the energy in a symbol is increased or the prior 
probability is changed. 


Some observations: 


e If Pria,| = Prla;] Vk,l, the decision rule is the minimum distance decision rule 


of the ML detector, and an estimate of the noise power is no longer required. 
e For large o”, the decision rule is heavily weighted by the prior probabilities. 
e The M decision regions are separated by (N — 1)-dimensional hyperplanes. 


e For the general case, with no restriction for equal energy signals or equal prior 
probabilities, there may be decisions which are not made no matter what the value 


of y. If either of the following conditions are true, ata, = af a, WI, k or Pr(ax) = 
k i 
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Figure 6.1: Decision regions for R? 


Pr(a;) #0 Vl,k, then there will always be a value of y that will lead to decision 


Qj. 


3. Probability of Error 

In this section the probability of field error for the MAP detector is derived, then 
a useful bound is developed that is more practical to calculate. An error occurs when the 
decision @ = a, is made when q;, 14k was sent. The probability of error is described 


by the following equation: 


Pr(error) = S- S° Pry € Rz|a;) Prlay| 


k Itk 


=O f rlylen) Prfou)ay 


(6.19) 


88 


where FR; is the region where decision k is made. The probability of error can also be 


written as 1 — Pr(correct): 


Pr(error) = 1 — Pr(correct) 


a S- Pr(y € Ry|@x) Pr(a) (6.20) 


=e d i. Pr(y|ax) Pr(ax) dy 


Both (6.19) and (6.20) involve multiple (N — 1) — fold integrals over decision regions, 
resulting in complicated integration limits. A simplification can be found by considering 
the pairwise union upper bound. 


To develop the upper bound, begin by rewriting the probability of error as 


Pr(error) = S- S- Prigi(y) > ge(y) |x| Pr(a@x) — all the joint terms (6.21) 
k I¢k 


where each Prigi(y) > gx(y)|ax| term is the pairwise probability of deciding @ = a, over 
Q@ = a, given a, was sent. The joint terms are the intersecting areas that get accounted 
for multiple times from only considering pairwise decisions. 

Calculation of Prigi(y) > ge(y)|ax] still involves an (N — 1) — fold integral, but 
the complexity of the decision region boundary is reduced to that of a single hyperplane 
defined by (6.16). The simplified decision region, along with the Gaussian form of p(y|a@), 
allows the transformation of the problem to a one-dimensional integral. 

We can now use the equation for the received vector y (6.3) in the decision boundary 
defined in (6.16) to obtain 








1 1 P 
+ 4Ple.-a)= fafa Jofao 
P 
Pr(ai) 





= —(a, — ay)" (ay — ay) — 207 In 


eles 


The noise vector z from Equation (6.3) is a vector of zero-mean i.i.d. Gaussian random 


variables with covariance o?I. The random variable 2(a; — a;)'z, is therefore zero-mean 
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with variance 07 = 407(a;, — @,)" (az — a)).' The pairwise probability of making decision 


Q@); given a, was sent is 


CCT eo fe ula aa 
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Where Rj, is the region for decision @ = a, when only two decisions are considered; the 





other decision being @ = a;. The term dj, is given by dj, = J (aK — a)" (ap — a). 
This result is used in (6.21) to obtain the union upper bound? for the probability 
of field error Pre(m Ap): 


Plane = es ax) > Q (3 ! 7 In Eee} (6.24) 


xk 





For high SNR the “joint terms”, neglected in developing this bound, become small and 
this bound becomes tight. 

A less complex upper bound may be found by considering the minimum distance, 
dmin, a8 a Worst case. Form the M-ary upper bound [8] by considering all of the sequences 


equally likely and recognizing that all of the distance terms dj, in Equation (6.24) will be 





‘Given a random vector with covariance matrix R, the vector y given by the linear transformation 
Az has autocorrelation function Ry = AR, A** [103]. 40?(a% — ax) (ax — ar) 
The union upper bound says Pr (J; A;) < 55; Pr(Ai) [11] 
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greater than or equal to dmin. The M-ary upper bound is then given by 





6.25 
as (6.25) 


where dmin is the minimum Euclidean distance between any of the field values. Notice that 
for this bound there are still only M possible field values transmitted rather than 2%, but 
the field values are assumed equally likely. Although this bound is simple to calculate, it 
is not very tight for large values of M. Hereafter, this bound (6.25) will not be used. 


4. Matrix Formulation 

For the purposes of simulation of field error rate and calculation of the upper bound, 
it is convenient to write the above equations in matrix notation. By writing the equations 
in a matrix format, some of the loops can be eliminated allowing the use of optimized 
matrix libraries. 

Consider the transmission of L fields. The received signal Equation (6.3) can be 


extended as a matrix of received fields as 
Y=B+4+Z, (6.26) 


where B is an N x L matrix of the L particular transmitted column vectors a, and N is 
the number of symbols in each vector a. Each column vector @ is drawn from the M < 2% 
values for the vector a with non-negligible a priori probabilities (p, > €). The matrix Z 
is an N x L matrix of the L noise column vectors z. 

The discriminant functions given in 6.11 can be calculated using matrix notation 
as follows: 


i] 
G = 7 In(po;,) — 5 On Aaa)! Of + A'Y), (6.27) 


where In(-) of a matrix is an element-by-element natural log. The vector 0; is a length j col- 
umn vector of ones. The vector p is an M x1 vector of the a priori probabilities for each pos- 


sible field Pr(@,) = py. The matrix Aga is an N x M matrix of column vectors, where the 


‘i 
i” column of Aga» is given by lai [OjaF[0] a,{lja*[1] --- a [N —-llaz[N — 1] . The ma- 


a 


trix A isan N x M matrix of the possible transmitted signals A = | Q, ... Ay-1 ' 
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Table 6.1: Example distributions for a 
field field values separation of 1 field values separation of 10 px In(px) 








Qo n-d n-50 0.011 -4.4 
a, n-4 n-40 0.003 -5.7 
Q2 n-3 n-30 0.005 -5.4 
az n-2 n-20 0.015 -4.2 
a, on-l n-10 0.109 -2.2 
Q, Uu n 0.519 -0.7 
QA nel n+10 0.291 -1.2 
Q7 n+2 n+20 0.031 -3.5 
Qg nd+3 n+30 0.004 -5.4 
Qg nu+4 n+40 0.003 -5.9 
Qi9 n+5 n+50 0.008 -4.8 








The resulting I x L matrix G can be used to determine the field-based MAP 
estimate for a. The maximum value in each column of G is associated with the estimate 


for a from the MAP decision rule. 


5. Results 

Consider the case of 11 likely field values shown in Table 6.1. Each field is 28 bits 
in length, (N = 28). With 28 symbols per field, there are 27° possible field values; for this 
example we consider M = 11 of them to be likely. In practice this distinction is made by 
considering the possible field values to be those whose a priori probability p;, > €, for some 


selected threshold «. By making this selection there is a residual error rate of 


je (6.28) 


kVpp<e 


The amount of tolerable residual error will ultimately determine the maximum practical 
value for €. 

Each field has a priori probability and values as shown in the Table 6.1. This 
section presents results for two cases. In the first case, the likely transmitted values are 
represented as consecutive numbers. In the second, the values are separated by steps of 
ten. In each case, the transmitted symbols are mapped to the field values using binary 
encoding. In each case, the Euclidean distance between field values depends on a nominal 


initial value n. 
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Figures 6.2b and 6.3b show the field error performance of the field-based MAP 
detector compared to the field error performance of a ML symbol-by-symbol detector for 
a field value separation of one and ten respectively. The performance of the two detectors 
is found by simulation. The upper bound is from Equation (6.24). This bound and the 
coherently detected antipodal signaling curves are analytic results. 

Figure 6.2b shows the field error performance for the field-based MAP detector 
versus the field error performance for a ML symbol-by-symbol detector. In this example 
the values are separated by a value of one, with distribution illustrated in Figure 6.2a. 
Notice in this case the field error rate of the field-based MAP detector is close to that of 
the bit error rate for the optimum detector for antipodal signaling. The field-based MAP 
detector performs significantly better compared to the field error rate performance of the 
symbol-by-symbol detector; providing approximately 2 — 6 dB of gain. 

Figure 6.3b shows the field error performance for the field-based MAP detector 
versus the field error performance for a ML symbol-by-symbol detector. In this example 
the values are separated by a value of ten with distribution illustrated in Figure 6.3a. 
Because of the increased separation between field values, the field error rate performance 
is superior to that of the bit error rate performance of the optimal detector for antipodal 
signals. This results in better field error rate performance than the corresponding symbol- 
by-symbol detector; providing approximately 6 — 8 dB of gain. 

The field error performance of the ML symbol-by-symbol detector is given by 


Pees 1a Pe), (6.29) 


where P,. is the probability of individual symbol error. The field error upper bound is the 
union upper bound (6.24). 

The field values used in the bounds and the performance curves are generated by 
randomly selecting a number n € {M/2,2" —M/2-—1}, and then generating the appropriate 
M = 11 field values according to Table 6.1. Assuming a symbol can take on one of two 
values, a field of length N can represent 2” different values. These 2% values are mapped 
to the integers from 0 to 2% — 1. It should be apparent that the distances dj, depends on 
the values a field can represent. 

Notice the significant difference in performance for the two cases. The reason for 
the difference in performance between the two is the pairwise distances for the data used 


in Figure 6.3b are greater than the pairwise distances for the data used in Figure 6.2b. 
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Figure 6.2: Field-Based MAP for 28 Bit Field Values Separated by 1 
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Figure 6.3: Field-Base MAP for 28 Bit Field Values Separated by 10 
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Consider the case N = 10, M = 4 and n = 20 and a field separation of ten. This will 
lead to the pairwise distances shown in Table 6.2. Now consider the same case with a 
field separation of 1. This leads to the pairwise distances shown in Table 6.3. Notice the 
difference in the minimum distance, d,,j;,, for the two cases. The field error performance 
of the receiver will increase as the Euclidean distance between the pairs increases. Assum- 
ing equally likely field values, the smallest pairwise distance dj, will have the greatest 


contribution to the probability of error. 


Table 6.2: Distances separation of 10 
pair Euclidean Distance 


digo 1.2649 
dio30 0.8944 
dio.4o 0.8944 
dyoo39 0.8944 
dyp.4g 1.2649 
doa 1.2649 








Table 6.3: Distances separation of 1 
pair Euclidean Distance 
dig29 1.0954 
dig21 0.8944 
digo. 0.8944 
dy9.21 0.6325 
dy9.22 0.6325 
da1,22 0.8944 








B. JOINT MAP DETECTOR (SYNCHRONOUS SIG- 
NALS) 


In the next chapter, Chapter VII, the Field-based MAP detector is extended to 
handle multiple asynchronous cochannel signals. Below the Field-based MAP detector is 
extended to the special case of synchronous cochannel reception. Synchronous reception is 
when the symbols from each source have identical start and stop times from the perspective 


of the receiver. 
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Consider two signals each with fields of length N, synchronized in time. As before, 
the k” field from the 7“ signal can take on one of M;,, values. The simultaneous syn- 
chronous reception of signal 0 and signal 1 with Mo, and MM, possible values respectively 
will, in general, lead to Mo, x Mi; combinations. Note the implicit assumption that the 
underlying waveforms used to generate the symbol sequences have possibly unequal signal 
powers, and a cross-correlation coefficient p, given by (6.30), with absolute value less than 
1, |p| < 1. Now the decision problem is deciding which of Mo, x My, possible field pairs 
was sent. The performance of the resulting detector will depend on the possible values of 
the fields for each signal. The pairwise distances between these field-pairs may be much 


smaller than the distances for the single case. 


1. Model 
Consider the synchronous joint reception model defined in Chapter II. Recall the 


cross-correlation coefficient between the two signaling waveforms defined as 


7 [ vo(twr(t) at 
pP = . 
VS voltualt) dty/f dn(twr(t) at 





(6.30) 





For notational convenience, let y’ be an interleaved vector of each matched filter output. 


For the two channel case: 


y’ = [volo] vsf0) voll) wl] --- wolV-a wfv-1} 621) 


In this section the prime is used to represent an interleaved vector. The vectors yp) and y, 
are the outputs of unit-energy filters ([ h(t)h*(t) dt = 1) matched to the first and second 


signal respectively, and sampled at the symbol period. The output vectors are 

Yo = Ao + pay + Zo, (6.32) 
and 

Y, = pA + A, + 21, (6.33) 


where the vectors ap and a are the transmitted values for signal 0 and signal 1 respectively, 
and the elements of z and z, are zero-mean i.i.d. Gaussian random variables each with 
variance o?. Note that the variance of the random variables contained in the vectors Z 


and z, are equal because they are both obtained from filtering the same AWGN process 
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with unit-energy filters. The cross-correlation of z) and 2, is given by E{zozj7} = o7 pl, 
where p is the cross-correlation coefficient. 

The covariance matrix of the noise component of the vector y’ defined in Equa- 
tion (6.31) is © given by 


bp 0 0 0 
p< BuO 0 0 
001 p 0 0 

S=o07|0 0 p 1 0 0 (6.34) 
0 da, wt WG 
O00 OOo «ch Spe 

The interleaved observation vector from (6.31) can be re-written as: 

y =a,+2 (6.35) 


where the vector aj, is an interleaved vector of one of the Mo,; x M1, possible transmitted 











pairs: 
9,5 [0] + pos, [0] 
pa,5[0] + a17(0] 
ao,j[1] + poral] 
Qa, — P,5 [1] aie Q1 1 1] (6.36) 
ao,j|N = 1] + PQA N-— 1] 
PQ, j [NV = 1] + Q1 1 N- 1] 
The vector aj, is built from one of the pairs of vectors described in 
{ (00, O10), (Q@o0,011), ***, (@0,0,@1,m,-1), (01,010), **°, (o,mMo-1; ousn-1)} 
(6.37) 
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The conditional PDF of the length 2.N vector y’ is then multivariate Gaussian given 
by 








Hal) = l ex —(y’ 7 ay) 7 (y! 7 a.) 
Plu lal) = Grama ( : ) 8) 


2: Discriminant Functions 
The discriminant function, as introduced in Section VI.A.2, for the conditional PDF 


given in Equation (6.38) is given by [102] 


(y! — Bae = 2) 





a(y’) =In(Pr{ai]) — ( 
(6.39) 





_ ny (YEE ty! = 2(aj) TE"! + (aj) TS *ee 
=in(Pr[ai) ~ ( : ) 


The second line of (6.39) follows from the identity S~1 = (X~1)*? from the complex 
conjugate symmetry of the covariance matrix. The term (y’ — a)’=~'(y’ — a’) is the 
square of the Mahalanobis distance [11] for y’ and aj. The terms y7h-!y’ are common 
to each discriminant; they may be removed leading to new a simplified discriminant: 
1 

ay’) = n(Prle’]) — 5(ai)"E tay + (ai)"D"y! (6.40) 
As in the single-channel case, the decision boundaries are still hyperplanes. Unlike the 
single-channel case, the hyperplane defining the decision boundary will not, in general, be 


perpendicular to the vector connecting a; and ai,. 


3. Probability of Error 

The development of the probability of error calculation closely follows that for the 
single-channel case developed in Section VI.A.3. The union bound is used in an identi- 
cal fashion. The only difference is in the calculation of the pairwise probability of error, 
Prigi(y’) > gx(y’)|a,|. To derive the probability of error, begin by writing out the descrip- 
tion of the decision boundary g;(y’) = gx(y’). 


I\Ty—-1,/ 
a) ay 


In(Pr{ai}) — { ; O,) Og 


+ (af)"D-1y! = In(Prlag)) — SS 8 4 (ayy! (6-41 
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then collect the y’ terms on the right hand side. 


(ay)7 Oa (a,) 7D, 
2 2 

Prlay}\  (a)*2 ta; (a) 7D tay, 

Prla’] ] | p) p) 


(aj)? Uy! — (ay)? Uy! = — In(Pr[ai]) + In(Pr[ay)) 4 








(al)"Ey — (ai)PEy’ =n ( 
(6.42) 


The decision boundary is then found by then assuming a, was sent and replacing the 


interleaved observation vector, y’, with the right hand side of Equation (6.35). 











! 2 Pria’ (a! )?S“ a! (a’ )To-1ay! 
(aj)"= ‘(arr + 2") = (a’,)7 = "(arr fe zt) —|n (Fe | l ; i ke ; : 
(aj, — a4)? E7t2' =In Prlaj)\) (ot Ear (aj) Eto 

Pr[aq| 2 , 5 





(6.43) 


Since z’ has covariance matrix X, the random variable (a) — aj,)’S~+z’ has variance 


of = (aj — a4,)"5(al — af)" 





3Given a random vector x with covariance matrix R, the vector y given by the linear transformation 
Ax has autocorrelation function Ry = AR,A*" [103]. 


(a — oi, PETE (af — a, PEN"? = (a - a )TEM (a - at.) 
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Following the same methods used in the single-channel case in Section VI.A.3 the 


probability of pairwise error is given by 


Prlaty’) > outu'lai = flutes) dy! 




















Ri, 
Pria;,| d?_ (1, k) 
=p ee Ty-l , = l k M\"% 
[le a’,) erie n Priav! 5 
ig (Fa | d2 (1k) 
P roas Ty 1 he Pr[Q; 2 
r Cc Q’,) Zz Ae 
Prl(QJ\ | dar (Lk) 6.44 
a (a — a’! Patz! | In (Fa + ( ) 
"| dar(Ls) S dur(1, Fi) 
PrlQi}\ , dig(Lk) 
=Q In (Tray) + 5 








Where dj,(1,k) is the squared Mahalanobis distance [11] between aj and ai, given by 
diy(l,k) = (xj — 0%) "E(x — a) (6.45) 


Notice if © = 07I then d3,(l,k) = d?,/o7, and Equation (6.44) reduces to Equation (6.23). 
Following the method of Section VI.A.3, the result in (6.21) can be used to obtain 


the union upper bound for the probability of field error for joint signals Pre jm ap): 


rlat]\ — du(l,& 
Preamp) 2 Pr(a's) ee. (az i In (Sa) | : ‘) (6.46) 


ixk 





4. Matrix Formulation 

As in the single-channel case, it is often useful to develop a matrix formulation 
of the equations to allow the use of vector math libraries for computation. The matrix 
formulation is developed for Equation (6.40), in a similar fashion to that of the single- 
channel case (6.27). 
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The received signal (6.31) can be extended as a matrix of received fields as 
Y= B’+2' (6.47) 


Where B’ is an 2N x L matrix of the L transmitted column vectors a’, and 2N is the 
number of symbols in each vector a’. Each column vector a’ is drawn from the Mj)M, < 
2? possible values for the vector a’. The matrix Z’ is an 2N x L matrix of the L noise 
column vectors 2’. 

The discriminant functions given in (6.40) can be calculated using matrix notation 
as follows 


1 
G = In(po;) — 5 (02 Aas-tar) OL H(A) sy) (6.48) 


The vector 0; is a length j column vector of ones. The vector p is an (MoM) x 1 vec- 
tor of the prior probabilities Pr(a@’;). The matrix Ags-1g* is a 2N x (MpM;) matrix of col- 
umn vectors, where the i” column of Ayy-1q« is given by diag{as=tal*}. The matrix A’ 
isan NxM matrix of the possible transmitted signals A’ = jax a’, QMo—1)(My—1) 
The resulting (/)M,) x L matrix G’ can be used to determine the field-based MAP 


estimate for a’. 


5. Results 

Results are presented below for the Joint field-based MAP detector with two syn- 
chronous signals. The values transmitted for each signal are determined using Table 6.1 in 
the same manner as for the single-channel case. 

Figures 6.4 and 6.5 show the field error performance of the joint field-based MAP 
detector compared to the field error performance of a Joint ML symbol-by-symbol detector 
for a field value separation of one and ten respectively. The performance of the two 
detectors was found by simulation; the upper bound and the coherently detected antipodal 
signal curves are analytic results. 

Figure 6.4 shows the field error performance for the joint field-based MAP detector 
compared to the field error performance for a Joint ML symbol-by-symbol detector. Note 
that a single-channel receiver for cochannel reception is not considered; the performance 
of a single channel receiver is unacceptably poor for highly correlated signals of similar 
received power levels. In this example the values are separated by a value of one with 
distribution illustrated in Figure 6.2a. Field error performance curves are generated for 


three different values of the cross-correlation coefficient, p. Notice the case of p = 0 
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gives identical performance to that of the single-channel case presented in Figure 6.2b. If 
p = 0, the underlying waveforms used to transmit the information are uncorrelated and 
the optimal joint receiver reduces to two single-channel receivers. 

Figure 6.5 shows the field error performance for the field-based MAP detector versus 
the field error performance for a Joint ML symbol-by-symbol detector. In this example the 
values are separated by a value of ten with distribution illustrated in Figure 6.3a. 


The field error performance of the Joint ML symbol-by-symbol detector is given by 
Pre =1—(1— Prue)” (6.49) 


Where P,, is the probability of individual symbol error from a ML detector. 

The field error upper bound is the union upper bound (6.46) developed in Sec- 
tion VI.B.3. The probability of bit error for antipodal is also included for reference. As in 
the case for the single-channel example, the performance will depend on the distribution 
of possible transmitted values for each channel. 

Notice the significant difference in performance for the two cases. The reason for 
the difference in performance between the two is the pairwise distances for the data used 


in Figure 6.5 are greater than the pairwise distances for the data used in Figure 6.4. 


C. COMPUTATIONAL COMPLEXITY 


Here, we briefly examine the computational complexity of a brute force implemen- 
tation of the MAP field detector with field length N. Calculation of the discriminants in 
Equation (6.11) involves 2N multiply and accumulate (MAC) operations. There are 2% 
discriminant functions resulting in a complexity of = 2’2N MACs per field. For a field 
length of 10 bits this amounts to approximately 20,000 MACs, for a field length of 15 this 
is on the order of 10° operations. There is an exponential dependence on field length that 


makes this problem unfeasible for fields consisting of many more than ten symbols. 
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Figure 6.4: Joint Field-Based MAP for 28 bit Field Values Separated by 1 
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Figure 6.5: Joint Field-Based MAP for 28 Bit Field Values Separated by 10 
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One obvious reduction can be made by setting field probabilities that are negligible 
to zero. If only M non-negligible field probabilities remain, the complexity is reduced to 


2MN. Even with M in the thousands this remains a tractable problem. 


D. CHAPTER SUMMARY 


In this chapter a field-based MAP detector is developed for a vector of symbols 
where the a priori probabilities are for a field of symbols rather than the individual symbols. 
A gain of 2 — 8 dB in field (or sequence) error rate performance is found when a priori 
knowledge is used compared to a symbol-by-symbol ML detector. The symbols are assumed 
independent, i.e., no ISI and no inter-symbol dependence introduced at the transmitter. 

The field-based MAP detector is then extended to multiple signals, where the fields 
are restricted to be synchronous. The synchronous joint field-based MAP detector is 
demonstrated to have performance gains of 3 — 18 dB when compared to the joint ML 
detector. These gains are very dependent on the a priori distributions. 

The next chapter develops the Joint Field-based BCJR algorithm for joint field- 
based MAP detection using field-based a priori information. This algorithm allows the 
computation of a field-by-field MAP estimate where asynchronous cochannel signals in- 
troduce a dependency between consecutive fields. It will be seen that in order to make a 
decision on a particular field, all of the received data must be examined, not just the data 


for the current field. 
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VII. JOINT FIELD-BASED MAP DETECTION 


In the first half of Chapter VI, the Field-Based MAP Detector is developed for 
an individual field of symbols. In the second half of that chapter, the Field-Based MAP 
Detector is extended to the Joint Field-Based MAP Detector for the special case of syn- 
chronous cochannel signals. For the special case of synchronous cochannel signals, the 
optimum MAP detector is developed by observing only the current field. This chapter 
extends the work from the second half of Chapter VI to the general case of asynchronous 
cochannel signals. For asynchronous signals, the optimal MAP decision can not be made 
by observing only the current field; the entire received vector must be used in order to 
make optimal MAP decisions. 

The development of the Joint Field-Based MAP detector begins with a model of 
the multiple signal sources in a form that is compatible with the BCJR algorithm. The 
BCJR (or MAP algorithm) is introduced in Chapter II and detailed in Appendix B. The 
BCJR algorithm provides a computationally efficient method to calculate the a posteriori 
probabilities of the state transitions of a Markov chain based on an entire received vector 
of decision statistics. The received signal is passed through a bank of matched filters, one 
for each waveform, to generate a vector of decision statistics. The matched filter bank is 
followed by a whitening filter which whitens the noise component of the decision statistics, 
and creates a causal relation among the decision statistics. The whitened decision statistics 
are then passed to a modified Field-Based BCJR Algorithm. The modifications to the 
BCJR algorithm involve the generation of a time-varying trellis, modified branch metrics, 
and a modified traversal method. Finally the FER performance results for several signal 


configurations are presented. 


A. SIMILARITY TO JOINT MLSD 


The signal model and trellis structure for the Joint Field-Based MAP Detector using 
BCJR is similar in many ways to the Joint MLSD [2] discussed in Chapter I. In both 
cases a set of matched filters generates decision statistics for each channel such that the 
decision statistics for channel 2 are, in general, correlated with the decision statistics for 
channel j. In the description of the Joint MLSD, based on [73] in Chapter II, a modified 


branch metric is used, thus eliminating the requirement for a whitening filter. The use 
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of a modified branch metric is known as the Ungerboeck approach [54,55]. For the Joint 
Field-based MAP detector it is not immediately obvious how to develop the equivalent 
modified branch metrics. Note that [104] does develop the Ungerboeck model for an ISI 
channel, but it is still not obvious how to extend this to the time-varying Markov source 
that results from cochannel reception. For this reason, a whitening filter is used, similar 
to the Forney approach [51,55]. The whitening filter has an added benefit of creating a 
causal discrete-time representation. The output of the whitening filter at sample n is a 
function of the current input at sample n, and N, — 1+ N, previous input samples, where 
N, is the number of signals and N, is the channel memory length, assumed here to be 
zero. As discussed in Chapter II, the computational complexity of traversing the trellis is 


exponential with the sample memory, i.e., exponential with (N, — 1+ N.). 


B. JOINT FIELD-BASED MAP AND BCJR 


As in the previous chapter, a posteriori probabilities are desired for each possible 
received field value. Selecting the field value with the maximum a posteriori probability 


minimizes the probability of field error: 
a, = arg max|P(axly’)], 
ak 


where ax represents a field value from the k“” signal and y’ is the vector of interleaved 
matched filter outputs. If the source model can be constructed in such a way that it can 
be written as a Markov chain, then the BCJR algorithm can be used to efficiently perform 


this calculation. 


C. MODEL 


Figure 7.1 illustrates the reception of two cochannel signals in AWGN. The symbol 
vectors @, and Qe, from signal 1 and 2 respectively, are input to pulse shaping filters 
Ui(t + 71) and we(t + 7). The output of the pulse shaping filters along with AWGN is 
received at the input to a bank of matched filters. Each filter, h, and ha, is matched to its 


respective pulse shaping filter w,,. The correlated outputs of the two matched filters are 
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Figure 7.1: Signal Model 


whitened with a whitening filter, and the resulting whitened decision statistics are passed 
to a decision algorithm, the Joint Field-Based MAP Algorithm. 

To extend this model to a GMSK signal, as used in AIS, one or more of the Laurent 
pulses (Appendix A) will be used for ~. The Laurent pulses have support greater than T,. 
This is a relatively straightforward extension but is left out to simplify the development. 
For an example see [68]. 

Consider the reception of N, signals in AWGN: 


Ns—1 


r(t)= S> sp(t) +20); (7.1) 


k=0 


where s;(t) is the complex baseband representation of a signal, and z(t) is a complex 
AWGN process with power spectral density No. Each signal s;,,(t) can be represented as a 


linear combination of amplitude modulated pulses 


sp(t) = » Axorg[i]e(t — iT — Te), (7.2) 


where the unit-energy signal pulse, ~;,(t), is only nonzero on the interval |0,7;], and 7; is 
the time offset for signal k, where by convention 0 < 7; < tT, < T;, Vj < k. The amplitude 
of the k'” signal is represented by Ag. 
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The signaling waveforms v, are, in general, not orthogonal to each other. Let p 


describe the inner-product of two synchronous (7; = 7%) signaling waveforms. 


p= f wyltugtnar (7.3) 


Let p;,, denote the cross-correlation between signaling waveforms corresponding to 
the correlation associated with symbol n of signal k and symbol n of signal 7, where by 


convention T; < Tp. 
pix =f vi(Ouilt — Arar. (7.4) 


This corresponds to cross-correlation between the current symbol index of each signal 
(regions a and d in Figure 7.2). 

Let p,,; denote the cross-correlation between signaling waveforms corresponding to 
the correlation associated with symbol n of signal 7 and symbol n — 1 of signal k, where 


by convention 7; < T,. 
Pk = / w;(t)p,(t + (LT, — Ar)) dt (7.5) 


This corresponds to the cross-correlation between the current and the past (or future) 
symbol index of each signal (regions c and 6 in Figure 7.2). 


The outputs of the set of matched filters at the sample time are 


mld] =f r(Ouglt - iT, — na, (7.6) 
where y;, are the unit-energy! impulse responses of the underlying signal waveforms used to 
generate each signal, s,. The matched filter outputs (yo, y1,---, YN,—1) can be represented 
as a linear combination of the individual responses from each signal: the desired response, 
plus the response of 2 x (N, — 1) symbols from the other channels. This is illustrated 
in Figure 7.2 for the case of three signals. In the example, the desired response is the 
second symbol from the second signal a1]. The response of the matched filter has a 
contribution from the previous symbol (a3[0]) and current symbol (a3[1]) from Signal 3, 
and a contribution from the next symbol (a[2]) and current symbol (a;[1]) from Signal 1. 


These contributions are represented by the regions labeled a, b, c, and d in Figure 7.2. 





'The pulse yx(t) is unit-energy if [°° we(t)vx(t) dt = 1 
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Figure 7.2: Three Asynchronous Cochannel Signals 


The contribution of each signal to the output of an individual matched filter is developed 
below. 
Begin by expanding Equation (7.6) using Equation (7.1) and (7.2) to obtain 


CO Ns 
ln] =f | SOY Asaslidvs(t - iD, — 1) +208) | Rte = 7, — ny) at 
—Oo} jar 4 








rit) (7.7) 
Ns = 
=50 > Aes] / w(t — iT, — 7))Wk(t — nT, — Te) dt 
gat. |G —oo 


a i’ w(t — nT, — T,) z(t) dt 


From the finite support of w the integral in the double summation of (7.6) can be simplified, 


So 


if t<n-1 or i>n+l 
s 1 i Shan 
/ w(t — tT, — 7;)Y,(t — nT, — Th) dt = 4 py, if i=n (7.8) 
a pee WE pA khasns 1] 
pes UO TRAC HS aH 1 
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Then Equation (7.8) can be used to simplify the double summation in Equation (7.7) to 


obtain 


+30 Avcaln] f " aby(t)OE(t — (te — 7) dt 


j<k a 
Pj,k 
Ts 
+Y Ajasln+ 1] f wby( E(t — Ts — (te —1)) at 
j<k ee ES Ee 
= ay (7.9) 
+ 2 Asagln — 1] f wh; (t)we(t + T, — (tT, — 7;)) dt 
j>k Ne - , 
Ts 
+) Ayosln] f wy (t)wk(t — (tT — 7;)) dt 
j>k 5 Fe 
Pj,k 


where z;[i] is a complex Gaussian random variable with variance 0? = No. The Gaussian 
random process formed by interleaving the noise components from the output of each 
matched filter is given by z[i] = E [¢] za[2] zs[é] ... zw, [é] pa and has autocorrelation 
matrices [2] 
oR], iff =i+1 
T oho od 
{ele = 47% Oh (7.10 
o?R'(lj, iff =i-1 


0 otherwise 


where Rii — j] is the (NV, x N,) matrix with elements n,m given by [2] 


Ramli —j] = ie Wn(t — iT, — Tr) 0", (t — jTs — Tm) at (10) 
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D. WHITENED MATCHED FILTER 


To easily calculate the branch metrics, and to reduce the number of symbols upon 
which a decision statistic depends, the decisions statistics are passed through a whitening 
filter. Take the N, length-N vectors of matched filter outputs and interleave them to form 
a vector of length-(V,N) matched filter outputs 


£ 
y= [vol] ml]... yx,1l0) voll] wll}... vol —1) wlN -1]...yn,alN - 4] 
(7.12) 
The idea of the whitening filter is to perform a linear operation on the interleaved vector 
y’ that removes the correlation between the noise components. 


The interleaved vector of matched filter outputs is written in matrix form as 
y’ = RAa’ + 2’, (7.13) 


where A is the diagonal matrix of signal amplitudes. The vector a’ is the interleaved 


vector of length N transmitted symbol vectors 


a’ = [ald] ay[0] ... ay,—1]0] oll] ai[1] ... ao[NM-—1] ai[N -1]...an,-1[N — 1] 7 


(7.14) 
and z’ is the vector of interleaved correlated noise components of the decision statistics y’, 
with Efz'z!™} =o°R. The matrix R is the global N,N x N,N correlation matrix for the 


noise components of the matched filter outputs 


Ro] (Ril)? 0 0 0 
Ri} Ro) (Rp)? oo .. 0 

R=|0 Rij Ro) (RI) ... 0 (7.15) 
0 0 0 0 R(I] R(0| 


where the elements of the matrices R[0] and R[1] are defined in Equation 7.11. 
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For illustration, a 4 x 4 correlation matrix for two signals is given below in Equa- 

tion (7.16) 
1 pi2 O 0 
p2 1 pa OO 
O pa 1 pre 
Oo, SO": pug? 


R= (7.16) 


where /j2 and fo; are the cross-correlation coefficients defined in Equation (7.4) and (7.5). 

Because the correlation matrix R is positive definite, the Cholesky decomposition 
of R can be used to find an upper triangular matrix U, such that R = UU [103, 105]. 
The noise components of y’ are then whitened by applying the transformation U~! to the 
interleaved vector of decision statistics defined in Equation (7.13). This transformation 


results in a new whitened set of decision statistics y’, [106] 


y, =Uy’ 
=U“ RAa!’ +U~*2z! 
=U 'UU" Aa! + U* 2! 
=U" Aa’ + 2/, 


(a7) 


/ 
where Z,,, 


is the whitened version of 2’ with €{z/,z’,"} = 07I, where I is the identity 
matrix. Since the matrix U~' is also an upper triangular matrix [107]; the whitening 
transform can be thought of as a noncausal filtering operation. This transformation not 
only whitens the noise components of the matched filter outputs, but also generates a new 
vector of sufficient statistics y’, whose signal components are a function of the current 
and past N, — 1 symbols. The causal transform results because the matrix U7 is lower 
triangular, so the resulting vector is a function of current and past values rather than 
current, past and future values. The vector will only depend on N, — 1 past values, rather 
than all past values, because U7 is banded?. This is in contrast to the vector y’ whose 
elements are also a function of the future N, — 1 symbols. 

This transformation is a time-varying finite impulse response (FIR) filter on the 


matched filter outputs y’. The transformation is time-varying because the transformation 





2Ann Xn matrix A with elements aij is banded if aj; =0 Vi <i—ky anda; =0 Vj >1+ ke for 
some 0 < ky, kg <n 
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state 





branch 


Figure 7.3: Nomenclature for Trellis 


matrix U~' is, in general, not Toeplitz. The outputs of the whitened matched filter, Ys 


may now be modeled as a Markov chain received through an AWGN channel. 


E. GENERATION OF THE TRELLIS 


The noiseless portion of y/, can be modeled as a time-varying Markov chain. The 
trellis is a tool used to represent the Markov chain as time progresses. The trellis consists 
of N + 1 stages, where N is the length of the whitened interleaved decision statistics, y/,. 
Each stage consists of a number of states representing the permutations of the memory of 
past transmitted symbols. The memory of each stage is large enough to account for the 
ISI or signal overlap, and also to account for the a priori probabilities for the fields. For a 
binary alphabet, each state in each stage has no more than two branches connecting it to a 
state in the next stage: one branch for each of the two possible symbol values. Each state 
in each stage will then have one or more branches connecting to a state from the previous 
stage. This is illustrated in Figure 7.3. 

The states in the trellis are defined in such a way that the Markov property of the 


state transitions is preserved 
Pr[A[n+1] = glA[n] = p] = Pr[Aln+]] = glAln] = pn, ANY] = Pri, ---,A[0] = po]. (7.18) 


For field widths greater than one symbol, this results in a time-varying trellis; this is 


illustrated below. 
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Table 7.1: Signal 0 Structure 








Field Number Width Values / Prior Pairs 

0 iL (0, 0.5 ) fbOas) 

1 3 (000,0.2) (101,06) (111, 0.2 ) 
2 1 (0, 0.5 ) OSs) 

3 1 (0, 0.5 ) fb 705<) 


Table 7.2: Signal 1 Structure 
Field Number Width Values / Prior Pairs 








0 1 (1,05) (0,05) 
1 2 (01,0.5) (00,05) 
2 1 (0,0.5) (1,05) 
3 D (00,0.5) (11,05) 


Consider two signals with parameters listed in Table 7.1 and Table 7.2. Each signal 
contains a various number of fields; each field with a various number of symbols. For each 
field there is a set of possible values and associated prior probabilities p,. The signals are 
received at a single receiver with their respective time offsets, 7 and 71, with 7 < 7}. 
After the whitening filter, outputs at time n is a function of the previous N, — 1 symbols. 
The time-varying trellis is generated in such a way that the Markov property (7.18) is 
maintained. This property is maintained if all the overlapping fields of the current symbol 


are stored in the current state. 


1. Determining the Stage Memory 

The formation of the trellis begins by determining the required number of past 
symbols represented by a stage (the memory of the stage). This involves looking back 
from the current symbol at all of the overlapping fields. The basic algorithm is outlined 
in Table 7.4. The algorithm starts with the interleaved vector of whitened matched filter 
outputs y’,. Each value in y’, is associated with a transition to a new stage in the trellis. 
Each value in yf, is also associated with a specific field of a specific signal. The first symbol 
included in the memory for a stage is the first symbol from the set of overlapping fields. 
The last symbol to be represented by a stage is the symbol associated with the current 
element of y’,. This results in the M symbols represented by this stage. Next, the number 
of unique elements represented by the M symbols must be determined; this is the number 


of states. 


116 


As an example of determining the memory of a stage, consider Figure 7.4 based on 
the signals listed in Tables 7.1 and 7.2. 


e At Stage 4, corresponding to the reception of symbol a,[1] on the Stage 3 to Stage 4 
transition, the current field (f;[1]) does not overlap the first symbol in Stage 3 mem- 
ory (a,[0]), but the next field (fo[1]) does. Therefore, the first element represented 
by this stage remains the same, and the states will represent combinations of the 


following symbols: a, [0], ao[1], ai [1]. 


e At Stage 5, corresponding to the reception of symbol ao[2] on the Stage 4 to Stage 5 
transition, the current field (fo[1]) overlaps the first symbol in Stage 4 memory 
(a;|0]). Therefore, the first element represented by the stage remains the same. The 


states will represent combinations of the following symbols: a,[0], ao[1], a1[1], ao[2}. 


e At Stage 6, corresponding to the reception of symbol a,[2] on the Stage 5 to Stage 6 
transition, field fo[1] is still not complete. Therefore, the states will represent com- 


binations of the following symbols: a,[0], ao[1], ai[1], ao[2], a1 [2]. 


e At Stage 7, corresponding to the reception of symbol ao[3] on the Stage 6 to Stage 7 
transition, the current field (fo[1]) does not continue, and the next field (f([2]) is 
not a continuation. Therefore, only ao[3] is required in memory; the states will 


represent combinations of the following symbol ao[3}. 


Table 7.3 lists the symbols comprising the memory of each stage based on the two 
signals listed in Tables 7.1 and 7.2. 


Daa Determining the Number of States per Stage 

As discussed in Chapter VI, some field values have negligible a priori probability. 
Ignoring field values with negligible probability leads to a reduction in the number of states 
required when constructing the trellis. An algorithm for determining the number of states 
is shown in Table 7.6. The algorithm starts by determining the number of elements as 
outlined in Table 7.4. Continuing with the example above, observe from Table 7.3 that the 
values stored in a state can span multiple fields. To find the number of states, generate all 
possible partial field values based on the symbols of a field represented by the memory of 
the current stage. Finally, find the unique words after this masking, and assign a state for 


each unique word. 


ee 
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Figure 7.4: Example Trellis Using Fields Described in Tables 7.1 and 7.2 


Table 7.3: Number Of Elements Per Stage 








Stage Number Elements 
1 Ag|0 

2 a,|0 

3 ay 0 ao(1] 

4 ay 0 ao(1] a,[1] 

5 ay 0 ao(1] a,[1] ao(2] 
6 a4[0] ao{1] ay[1]  a@o[2] ay [2] 
7 ag 3 

8 at 3 

9 ao 4 

10 ao[4] a, [4] 

11 ao[4] ai[4]  ao[5] 

12 ay,|o 
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Table 7.4: Determining the number of elements represented by a stage 
for all Elements of y/, do 
CurrentFirstElement < Current Index 
OtherFirstElement < Current Index 
if Current Field Continuing then 
CurrentFirstElement <= Beginning of current field —1 
end if 
if Other Field Continuing then 
OtherFirstElement < Beginning of other field —1 
end if 
FirstElement <= min(OtherFirstElement,CurrentFirstElement) 
end for 


Table 7.5: Permutations of Field Values for Stage 4 









































Al] Al] fold] 
1 00 ~—000 
0 00 000 
1 01 000 
0 O01 000 
1 00 101 
0 oO 101 
D- “0, @0n 
0 ol 101 
Ti OO © S00 
0 oO 111 
t On 41d 
0 ol ttl 





For example, consider Stage 4 of Figure 7.4 representing the following symbols 
a[0}, ao{1], a: [1], this stage contains symbols from fields zero and one from Signal 1 (f;[0], 
fi{l]). Stage 4 also contains symbols from Field 1 of Signal 0 (fo[1]). The possible per- 
mutations of these three fields with the symbols a;|0], ao[1], a1{1] highlighted, are shown in 
Table 7.5 Although there are 12 permutations of field values, Stage 4 only requires four 
states; the other permutations of the fields are for symbols that are not contained the 


Stage 4 memory. 
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Table 7.6: Determining the number of required states for each stage 


for all Stages do 
Determine number of elements for this stage 
Determine the number of unique groups in this set of elements 
for all Unique Element of Set do 
Assign State Element of Set 
end for 
end for 


3. Making the Branch Connections 

Once the values represented by each state are determined, the next step is to make 
the branch connections. Table 7.7 summarizes the algorithm for connecting the states from 
one stage to the next stage. At each stage this involves a loop through all of the states at 
the next stage for each state at the current stage. Working from left to right across the 
trellis, at Stage n State p there are at most K branches leaving each state, where K is 
the size of the alphabet for the symbols. A particular state represents a permutation of 
the received symbols represented by this stage. A branch exists from State p in Stage n 
to State g in Stage n + 1 if the mutual symbols in both states are all the same value. If 
a particular value for the next symbol does not lead to a field value with a probability 
of occurring greater than ¢, (py > €), then the state transition will not be drawn. The 
resulting trellis with branch connections is shown in Figure 7.5. Consider the branches 
connecting Stage 8 to Stage 9 in Figure 7.5. There is only one branch leaving State 0. This 
is to be interpreted as follows: with the information in State 0 at Stage 8 there is only one 


possible value for the next symbol. 


F. BRANCH METRIC CALCULATION 


The branch metric calculation for the field-based BCJR algorithm requires a slight 
modification from the standard method outlined in Appendix B. The term ¥y represents 


the probability of transition from State p at Stage n to State q at Stage n + 1 with the 
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Figure 7.5: Resulting Joint MAP Trellis 


Table 7.7: Determining the Branch Connections 


for all stage do 
for all cstate in current stage do 
shift current values over one bit 
mask off overlapping elements 
for all nstate in next stage do 
if cstate in current stage contains same elements nstate then 
add branch connection 
end if 
end for 
end for 
end for 
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observation y,,. The term y can be factored as follows: 


Yn(P, g) =Pr(Aln + 1] = 9g, ynlAln] = p) 
=Pr(yplAln + 1] = 9g, An] = p) Pr(Afn + 1] = GlAln] = p) Gly) 


a priori probabilities 


where A|[n] is the state at Stage n. The second factor is the a priori probability Pr(a, = a). 
The first factor can be rewritten as Pr(yp|Sp) where spq is the output of the channel 
uniquely associated with a state transition from p to g. Because the noise components 


have been whitened, the conditional probability can be written as 


1 —llyn = Spall? 
Pr(yn|Spg) = (naa (“ape (7.20) 


where o? is the variance of the noise term from Equation 7.17 and ||(-)||? is the L2-norm. 

The modification for field-based BCJR is in the second factor, the a priori prob- 
ability Pr(a, = a). The algorithm is provided a priori field probabilities; a specific field 
probability is represented by the transition from stage to stage on the last bit of a field. If 
the state transition does not represent the last symbol in a field, then this probability is 
set to one. The a priori probability for the branch metric could arbitrarily be set to any 
constant value for each transition; the effect is removed during normalization. If the state 
transition does represent the last symbol in a field, then use the a priori probability for 


that field value [108]. 


G. TRELLIS TRAVERSAL 


Trellis traversal is identical to the standard BCJR algorithm detailed in Appendix B. 
For the forward pass, the values for a at Stage n and State g are calculated by the following 


recursive equation: 


Q-1 
Anti(g) =), Pr(Aln] =, Yen) Praln + 1] = @, wlAln] = p) 
—_—oO OO 
p=0 
an (p) Yn (P-4) (7.21) 
Q-1 
p=0 


de 


where Q is the number of states at Stage n. 
In a similar manner for the backwards pass, the Gs at Stage n and State p are given 


by the following recursive equation: 


Q-1 
Bn(p) = S— Pr(A[n + 1] = 9, ynlAln] = p) Pr(ysnlAln + 1] = 9) 
‘n(P.4) Bn+1(4) (7.22) 
Q-1 
=) nlp, ¢)Bn+41(q) 
q=0 


where @ is the number of states in the next state at time n+ 1. 


After calculating the as and the (s, the a posteriori probabilities are then given by 


Pr(a = az|y) = piess On(P)Yn(P, Dn+1(G)/ Pr(y) 
“ DareA La(p,q)eSe On(P)In(P, 2) Pn+1(@)/ Pr(y) 





(7.23) 


= Da (pa) eSe On(P)Y¥n(P; 1) On+1(q) 
ere wa Hes: An(P)Yn(P; 7) Bn+1(Q) 





where S, is the set of transitions associated with the specific input a, and a, € A where 
A is a finite alphabet. 

The effort to preserve the Markov property is only to keep track of the a priori 
field probabilities. For example, the received values in y’, used in the metric for the 
values for y (branch metrics) only depend on the last NV, — 1 symbols placed into the state 
and the current symbol. This is the memory from the channel, or equivalently from the 
asynchronous overlapping signals. The additional state memory beyond N, — 1 allows for 
book keeping of the a priori probabilities. A consequence of this is that many of the branch 


metrics (y) are identical; they are not be a function of the additional stage memory. 


H. RESULTS 
This section presents the FER performance results for the Joint Field-Based MAP 


Detector. Performance results of the detector are presented for two signals with a simple 


structure leading to the Joint MLSD developed in Chapter II. Results are then presented 
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for several other more complex signal configurations. Table 7.8 summarizes the cases 


investigated in this section. 


Table 7.8: Table Summarizing Test Cases 
Case Description 








A Simple case, equally likely independent symbols, each field is length one 
D 3-symbol fields with two equally likely values 0 and 1 
F 3-symbol fields with general distributions for field values 
B various field widths, various field value distributions 
1. Case A: 


In this first example the performance of the Joint Field-Based MAP Detector is 
demonstrated for the structure of the signals such that it reduces to a joint MLSD. 

The structure of the two signals are shown in Table 7.9 and Table 7.10. Each 
signal consists of single-symbol fields which take on one of two values with equal a priori 
probability. Notice that the structure of the trellis in Figure 7.6 is identical to the structure 
of the joint MLSD introduced in Chapter II. This simple structure results for the two signal 
case when the fields are only one symbol wide. 

a. Signal Offset Ar = OT, 

First, it is necessary to point out that here Ar = O7;, and there is no ISI 
introduced by the transmitter or channel, therefore there is no need for the trellis approach. 
The optimum receiver can be developed as in Chapter VI. The example of Art = OT, is 
used to demonstrate that for this case the algorithm performs equivalently to the method 
outlined in Chapter VI. In this first example FER curves are plotted (Figure 7.7) for 
field three of single one with various values for the cross-correlation coefficient p defined in 
Equation (7.3). 


Table 7.9: Case A: Signal 0 Structure 
Field Number Width Values / Prior Pairs 








0 1 (0,05) (1,05) 
1 1 (1,05) (0,05) 
2 1 (1,05) (0,05) 
3 1 (1,05) (0,05) 
4 1 (1,05) (0,05) 
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Table 7.10: Case A: Signal 1 Structure 
Field Number Width Values / Prior Pairs 








0 1 (1,05) (0,05) 
1 1 (1,05) (0,05) 
p i (1,05) (0,05) 
3 i (1,05) (0,05) 
4 i (1,05) (0,05) 


Table 7.11: Case A: Number Of Elements Per Stage 
Stage Number Elements 








ao 
ay 
ao 
ay 
ao 
ay 
ao 
ay 
ao 
0 ay 





rE ew WM NYOrF HY COO 





F OCON nok wd 


Notice that with p = 1 the FER is approximately 0.5 (actually it approaches 
0.25 for high SNR, of the four possible transmitted pairs only (—1,+1) and (+1,—1) lead to 
an ambiguous result that will likely result in an error). Notice that the BER is identical to 
the results presented in Chapter II for Joint Detection. For a cross-correlation coefficient 
of p = 0, the performance is identical to single-channel binary antipodal signaling. For 
small cross-correlation coefficient, the FER performance is very close to binary antipodal 


signaling. 













Figure 7.6: Case A: Trellis 
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Figure 7.7: Case A: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, synchronous case At = OT, 


b. Signal Offset At = 0.5T, 

With an offset between the two signals reasonable FER performance can be 
achieved even with highly correlated signals. This is similar to the example outlined in 
Section I.E. Although the two fundamental pulse shapes have a cross-correlation coefficient 
of p = 1, the time offset leads to pyg = po, = 0.5. Because of this, the signals can be 
resolved. Figure 7.8 illustrates this for the third field of signal one; even for a cross- 
correlation coefficient of 1.0 the FER performance is acceptable. 

There is a variation in performance based on the distances from the start 
of the packet; in order to make the decision on symbol n, all symbols received prior to n 
must not be in error (or there must be an even number of errors). Figure 7.9 presents FER 
performance for various field numbers; larger field numbers correspond to greater distances 
from the start of a packet. Notice that the error rate is higher for fields that are further 
away from the beginning of the packet. 
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Figure 7.8: Case A: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, asynchronous case At = 0.57, 
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Figure 7.9: Joint Field-Based MAP Performance Curves Variation Based on Distance from 
Start of Packet Correlation Coefficient = 1, asynchronous case At = 0.57; 
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2. Case D: 

In this example, the FER performance of the Joint Field-Based MAP Detector for 
signals with a fixed field width of three symbols is examined. Each field takes on one of 
two possible values, 000 or 001. We might expect the performance in this example to be 
very similar to that of the previous section,(VII.H.1); as illustrated below, this is not the 
case. The structure of the signals is shown in Tables 7.12 and 7.13. The structure of the 
resulting trellis is illustrated in Figure 7.10. Notice the simple state transitions that result 


from the signal structure. 


Table 7.12: Case D: Signal 0 








Field Number Width Values / Prior Pairs 
0 3 (001, 0.5) (000, 0.5 ) 
1 3 (001, 0.5) (000, 0.5 ) 
2 3 (001, 0.5) (000, 0.5 ) 
3 3 oe 0.5) (000, 0.5 ) 
4 1 COs) 


Table 7.13: Case D: Signal 1 








Field Number Width Values / Prior Pairs 
0 3 (001, 0.5) (000, 0.5 ) 
1 3 (001, 0.5) (000, 0.5 ) 
2 3 (001, 0.5) (000, 0.5 ) 
3 3 (001, 0.5) (000, 0.5 ) 
4 1 (0,1) 





10 11 #12 #13 #14 «15 «#16 «17 ~«218«19 «20 «#21 «22 «23 «24 «225 «(26 
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Figure 7.10: Case D: Trellis 
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Table 7.14: Case D: Number Of Elements Per Stage 





























Stage Number Elements 
1 ao 0 
2, ao 0 ay 0 
3 ao 0 ay 0 ao(1] 
4 ao(0] ay [0] aol] aif 
5 ag[O]} = ay[O]  ao[ 1] ai [1] a [2] 
6 ay|2 
7 ay,|2 ao|3 
8 ay|2 ag|3} ay }3 
9 ay[2] — ao[3] a1 [3] a4] 
10 a4[2] — ao[3] ay [3] ao[4] ay [4 
11 ay 2 ao 3 ay, 3 ao|4] a, [4] a5] 
12 ay|o 
13 ay,|9 ao|6 
14 ay,|o ao|6} a4 |6 
15 ai[5| do [6] a, [6] a7] 
16 at 1) ao 6 ay 6 ao|7| a,|7| 
17 ay[5] a6] ax[6]  aol7] aa [7] ao [8] 
18 ay 8 
19 a4|8 ao|9 
20 ay[8] — ao[9] a1 [9] 
21 ay[8] — ao[9] a1 [9] ao [10] 
22 ay 8 ao 9 ay [9] ao|10} ay 10] 
23 ay 8 ao 9 ay [9] ao[10] 1{10] ao[11] 
24 a,(11] 
25 ao(12] 
26 a(12] 
a. Signal Offset Ar = OT, 


Figure 7.11 presents the FER performance of the Joint Field-Based MAP De- 
tector for synchronous signals and various cross-correlation coefficients. The performance 
here is identical to the results for Case A presented in Section VII.H.1. 

b. Signal Offset At = 0.5T, 

Figure 7.12 presents the FER performance for Joint Field-Based MAP De- 
tector for the asynchronous case for various cross-correlation coefficients. The performance 


for this signal structure is superior to the results presented for Case A in Section VII.H.1. 
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Figure 7.11: Case D: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, synchronous case At = OT;, fixed width fields 


The performance is only a fraction of one dB from the antipodal signaling case. The reason 
the detector performs so well, even with p = 1, is that each of the three bit fields only has 
one unknown bit. With an offset between the two signals, Av > 0, the unknown bits in 
each field partially overlap the known bits in the other signal. 

In Figure 7.13 the FER performance for the Joint Field-Based MAP Detector 
are presented for the asynchronous case for various Fields. The dependence of error rate on 
distance from the beginning of a packet, in contrast with Case A in the previous section, is 
not seen here. For Case A there is a requirement that the previous n fields not be in error; 


in this case, because of the known symbols within the fields, there is no such requirement. 
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Figure 7.12: Case D: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, asynchronous case At = 0.57, 
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Figure 7.13: Case D: Joint Field-Based MAP Performance Curves for Various Field Num- 
bers Correlation Coefficient p = 1, asynchronous case At = 0.57, 
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3. Case F: 

In this next example the results for a detector with fixed width fields but variable 
distributions are presented. The signals, with structure illustrated in Table 7.15 and 7.16, 
have field widths identical to Case D presented in Section VII.H.2. As illustrated in 
Figure 7.14, this results in a much larger trellis. Figure 7.15 shows the same trellis enlarged 


and cropped for legibility. 


Table 7.15: Case F: Signal 0 








Field Number Width Values / Prior Pairs 

0 3 (001, 0.2) (000,05) (101,03) 

1 3 (000, 0.125) (001, 0.125) (010, 0.125 ) 
(100, 0.125.) (101, 0.125")> -(.110,.0.125. } 

2 3 ( 000, 0.01 ) ( 001, 0.01 ) ( 010, 0.01 ) 
(100,0.1) (101,06) (110,01) 

3 3 (000,0.5) (111,05) 


Table 7.16: Case F: Signal 1 








Field Number Width Values / Prior Pairs 
0 3 ( 001, 0.5 ) ( 000, 0.5 ) 
1 3 ( 001, 0.2 ) ( 000, 0.2 ) (OF1,0;2 ) 
2 3 ( 000, 0.125 ) (001, 0.125) (010, 0.125 ) 
C100.0.125 J (COL 0125.) «CTO 0.125.) 
3 3 ( 001, 0.5 ) ( 000, 0.5 ) 
a. Signal Offset Ar = OT, 


Figure 7.16 presents the FER performance for Signal 0 Field 3 of the Joint 
Field-Based MAP Detector for synchronous signals and various cross-correlation coeffi- 
cients. Note from Tables 7.15 and 7.16 that the two fields, Signal 0 Field 3 and Signal 1 
Field 3, are further in Euclidean distance than is the case for the signals presented earlier. 
This increase in distance reduces the probability of field error. 

b. Signal Offset At = 0.5T,; 

Figure 7.17 presents the FER performance for Signal 0 Field 3 of the Joint 
Field-Based MAP Detector for asynchronous signals with a separation of Ar = 0.57;,, and 


various cross-correlation coefficients. The performance here is slightly better than that 
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Figure 7.14: Case F: Full Trellis 
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Table 7.17: Case F: Number Of Elements Per Stage 





























Stage Number Elements 

1 ao 0 

2, ao 0 ay 0 

3 ao 0 ay 0 ao(1] 

4 ao(0] ay [0] aol] aif 

5 ag[O]} = ay[O]  ao[ 1] ai [1] a [2] 

6 ay|2 

7 ay,|2 ao|3 

8 ay|2 ag|3} ay }3 

9 ay[2] — ao[3] a1 [3] a4] 

10 a4[2] — ao[3] ay [3] ao[4] ay [4 

11 ay 2 ao 3 ay, 3 ao|4] a, [4] a5] 
12 ay|o 

13 ay,|9 ao|6 

14 ay,|o ao|6} a4 |6 

15 ai[5| do [6] a, [6] a7] 

16 at 1) ao 6 ay 6 ao|7| a,|7| 

17 ay[5] a6] ax[6]  aol7] aa [7] ao [8] 
18 ay 8 

19 a4|8 ao|9 

20 ay[8] — ao[9] a1 [9] 

21 ay[8] — ao[9] a1 [9] ao [10] 

22 ay 8 ao 9 ay [9] ao|10} ay 10] 

23 ay 8 ao 9 ay[9] ao[10] 1{10] ao[11] 
24 a,(11] 


of the synchronous case presented above. For this case the FER performance for each 
field will depend on the a priori distributions as outlined in Table 7.15 and 7.16. This 
variation is illustrated in Figure 7.18 for the synchronous case and in Figure 7.19 for the 


asynchronous case. 
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Figure 7.16: Case F: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, synchronous case At = OT,, fixed width fields 
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Figure 7.17: Case F: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, asynchronous case At = 0.57, 
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Figure 7.18: Case F: Joint Field-Based MAP Performance Curves for Various Field Num- 
bers Correlation Coefficient p = 1, synchronous case At = OT, 
































































10! - 7 
10° E : 
107 E 4 
10-? E : 
10-3 E Antipodal @ | 
ae ia chan 0, fieldO © +4 
10 : chan 1,field0 O F 
107° ‘ anl,fiedd1 + 4 
fool chan 0, field 3 : 
10°" chan 1, field3 @ + 
1078 £ | | | | | 1 
-5 0 5 10 15 20 25 


E»/No (dB) 


Figure 7.19: Case F: Joint Field-Based MAP Performance Curves for Various Field Num- 
bers Correlation Coefficient p = 1, asynchronous case At = 0.57, 
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A. Case B: 

In this final example the performance of the Joint Field-Based MAP Detector for 
signals with variable field widths is presented. 

The structure of the two signals are shown in Table 7.18 and Table 7.19. Each 
signal consists of fields of various widths which take on various distributions. Notice that 
the structure of the trellis in Figure 7.20 changes for the different field widths and prior 


distributions. 


Table 7.18: Case B: Signal 0 
Field Number Width Values / Prior Pairs 








0 1 (0,05) (1,05) 

1 1 (0,05) (1,05) 

2 1 (1,19 

3 3 (000, 0.2) (101,06) (111, 0.2) 
4 1 (0,05) (1,05) 

5 1 (0,05) (1,05) 

6 1 (0,05) (1,05) 


Table 7.19: Case B: Signal 1 
Field Number Width Values / Prior Pairs 








0 1 P08 Je Oe 05.) 
1 1 0S" .H02057) 
2 2 (01,0.5) (00, 0.5 ) 
3 1 (OOo). < C102) 
4 2 (00,0.5) (11, 0.5 ) 
5 1 CO0ibey bs O58.) 
6 1 COS0s ye = beOar) 
a. Signal Offset At = OT, 


In this first example, FER curves are plotted for one of the fields for the 
synchronous case and various values for the cross-correlation coefficient p (7.3). 
Figure 7.21 illustrates the FER performance for Chan 0 Field 3. The per- 


formance here is slightly better than some of the previous cases. 
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Table 7.20: Case B: Number Of Elements Per Stage 
Stage Number Elements 
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b. Signal Offset At = 0.5T; 

As in Case A, with an offset between the two signals we can achieve reason- 
able FER performance even in the case of highly correlated signals; Figure 7.22 illustrates 
this. Also, notice in Figure 7.23 there is not a large variation in performance based on field 


number. 
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Figure 7.20: Case B: Trellis 
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Figure 7.21: Case B: Joint Field-Based MAP Performance Curves for Various Values of 
Correlation Coefficient p, synchronous case At = OT, 
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Figure 7.22: Joint Field-Based MAP Performance Curves for Various Values of Correlation 
Coefficient , asynchronous case At = 0.5 
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Figure 7.23: Joint Field-Based MAP Performance Curves Variation Based on Distance 
from Start of Packet Correlation Coefficient p = 1, asynchronous case Ar = 0.57, 
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I. CHAPTER SUMMARY 


The Joint Field-Based MAP Detector is introduced in this chapter. This detector 
optimally incorporates a priori knowledge of variable width field values to aid in joint 
detection. The Joint Field-Based MAP detector can also be applied to situations where 
knowledge of one of the cochannel signals is very limited. If there is no knowledge about 
the field structure of one of the signals the knowledge of the other signals can still be 
exploited to aid in recovery. The signal model is developed before and after the whitened 
matched filter. We then developed the form of the trellis required to use the BCJR algo- 
rithm to generate the MAP estimates. Finally results are presented for nominal scenarios 
investigating how the performance varies as the structure of the signals change. These re- 
sults demonstrate the potential performance improvement that is possible by incorporating 


a priort information into a joint detection algorithm. 
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VIII. CONCLUSIONS AND FUTURE WORK 


The goal of this work is to investigate methods for efficient recovery of the in- 
formation content from low-entropy narrowband cochannel signals. There is significant 
prior work in joint detection of narrowband signals [2], but little work has focused on low- 
entropy signals. In this dissertation, optimal (minimum probability of field error) methods 
are developed to achieve this goal. A summary of the contributions from this work and 


suggestions for continued research are presented here. 


A. SUMMARY 


In Chapter II, the fundamentals of multi-user detection required for describing 
joint detection are described. A detailed background of the signaling waveform used in 
AIS (GMSK) is provided. A whitened discrete model of the signal is developed using 
the Laurent representation of GMSK. The MAP algorithm is then introduced as an effi- 
cient method for calculating a posteriori symbol probabilities. A literature review of joint 
detection methods is presented, focusing on work by Verdt [2]. 

In Chapter III, a complete overview of the AIS is presented. The signaling waveform 
and the packet structure of the AIS are presented. A link budget is presented demonstrat- 
ing that for nominal conditions, a receiver collecting from low earth orbit (approximately 
800 km) will not be noise-limited. An example of non-coherent AIS demodulation is pre- 
sented showing each step required to recover the transmitted data. The low-entropy prop- 
erty of the AIS message is demonstrated with an example. A 168-bit AIS packet is found 
to typically contain approximately 20 bits of new information. Finally the SOTDMA algo- 
rithm is investigated, and the predictability of subsequent AIS transmissions is presented. 

Chapter IV introduces methods of prediction focusing on AIS field values. Three 
different candidate predictors for the AIS fields are presented, and tuning parameters for 
the predictors are heuristically developed. Further, a novel predictor for the AIS commu- 
nications state field is presented, and a demonstration of the performance of a Kalman 
filter on AIS latitude and longitude fields. 

In Chapter V, MLSD is investigated for a bit-stuffed data source. The bit-stuffing 
process is modeled in a form compatible with a trellis representation. With this repre- 


sentation, the Viterbi algorithm may be used to find efficiently the most likely sequence. 
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The bit-stuffing process is found to be a catastrophic encoding process (a single bit error 
leads to a large number of decoder errors). Approximate BER bounds for MLSD of a bit- 
stuffed data source are developed. The MLSD is shown to have superior BER performance 
compared to symbol-by-symbol decisions (a gain of 0.2 — 1.0 dB for the unstuffed data). 

The Field-Based MAP detector is developed in Chapter VI. The performance of this 
detector is demonstrated using representative data. The detector is shown to outperform 
the symbol-by-symbol ML detector for some example a priori distributions. The detector is 
extended to multiple synchronous signals, and the performance of the resulting detector is 
presented. The performance gains from Field-Based MAP detection are highly dependent 
on the level of a priori information available and, in the joint case, on the correlation 
coefficient. The examples in Chapter VI demonstrate a performance improvement of 3 to 
18 dB of gain when using joint Field-Based MAP compared to joint symbol-by-symbol 
detection. 

Finally, Chapter VII develops the Joint Field-based MAP detector for asynchronous 
signals. The algorithms provided in this chapter illustrate the construction of a trellis 
compatible with the BCJR algorithm ultimately allowing for the calculation of a posteriori 
probabilities for transmitted field values. An efficient implementation of the Joint Field- 
Based MAP detector is developed and field error rate performance is investigated. The 
performance gains here are identical to that of the previous chapter for the synchronous case 
(Ar = 0). For asynchronous signals the results from Chapter VI are not applicable, and the 
detector from Chapter VII must be used. For the asynchronous case there are additional 
FER performance gains of up to 10 dB, depending heavily on the cross-correlation between 


the two signals. 


B. CONCLUSIONS 


Figure 8.1 gives a “big picture” view of the objective in this work as applied to 
remote reception of the AIS. The algorithm begins with some knowledge about existing 
ships — this could be very specific information such as a previous position report or very 
general, such as “Ship A left Port B at 1300 yesterday.” At the receiver, two or more 
signals are received overlapping in frequency and time. The signal parameters for each of 
the overlapping signals are estimated. This reception is passed through multiple matched 


filters, and each filter is matched to a specific pulse with estimated parameters. The 
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Figure 8.1: “Big Picture” 


resulting decision statistics are then passed to the Joint Field-Based MAP detector, along 
with a priori information derived from previous receptions or other sources. The output 
of the decisions from the joint detector are then passed to prediction algorithms for future 
receptions. This dissertation develops the important pieces for achieving this goal: the 
joint field-based MAP detector, the prediction method for the AIS fields, and the MLSD 
for bit-stuffing. 

The method depicted in Figure 8.1 is relevant even without a priori information 
for each of the simultaneous users. For example, suppose from previous transmissions it is 
known that Ship A will transmit on time slot 7. Based upon the reception at time slot 7, 
it is determined that there is an additional unknown ship, Ship X, occupying the same 
slot. With no knowledge of what Ship X is transmitting the field structure is simply set to 
equally likely symbols. The algorithm described in Chapter VII is now applied as before. 
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Many ideas from this work may be incorporated into receivers immediately. The 
work on MLSD in Chapter V is applicable whenever the source of the data has undergone 
bit-stuffing. If a receiver is already using the Viterbi algorithm for detection, the ideas 
from Chapter V can be incorporated by modifying the existing trellis. If the receiver is 
currently making symbol-by-symbol decisions, the Viterbi algorithm outlined in Chapter V 
may be implemented directly. The work in Chapter V demonstrates that a gain in BER 
performance of approximately 0.2 — 1.0 dB can be expected by using these techniques. 

The work in Chapter VI can be readily extended to any problem where a priori in- 
formation exists for sets of symbols. This is applicable to both single- and multiple-channel 
receivers. Both the gain in performance and the additional computational complexity will 
be very dependent on the particular situation. For example, if it is known that a set of 
symbols can take on very few values, then the computational complexity of implementing 
the field-based MAP detector is relatively low. If the set of values for a field of symbols 
are all separated by something greater than the minimum distance, then the performance 
improvement will likely be significant. If the set of possible values is very large, and the 
distance between pairs of values is approximately the minimum distance, then the small 
performance improvement may not justify the large computational burden of the field- 
based MAP detector. 

Chapter VII presents a culmination of the ideas from previous chapters. A method 
for calculating a posteriori probabilities for a transmitted field is developed. In order to 
implement this method of signal recovery the trellis must be set up according to the known 
structure and a priori information of the signals. The structure of the trellis will vary as 
the information available for the signals changes. For the general case of asynchronous 
reception the method illustrated in Chapter VII generates a posteriori field probabilities 
leading to the optimal (minimum probability of field error) detector. Gains of up to 10 dB 
are seen compared to the synchronous joint detector presented in Chapter VI. 

The theoretical work developed in Chapters VI and VII can be put to practical 


use as a bound for system performance. The detectors developed are optimal given the 
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stated constraints. The FER curves from these chapters can be used as an upper bound 


on performance. 


C. FUTURE WORK 


In performing this research, many areas deserving additional investigation have 
been identified. Most of the future work involves solving practical problems. The recom- 
mendations can be roughly broken into three categories: FER performance improvement, 
parameter estimation, and computation reduction. 

The work has been based on the first term of the Laurent decomposition. The trellis 
could be extended to use all the pulses, increasing performance. For simplicity, it is assumed 
there is no ISI between pulses; this is not true; GMSK introduces ISI. The performance 
demonstrated in this dissertation is based on a simplified channel. It is a straightforward 
extension to account for the ISI by adding states to the trellis. Furthermore, the MLSD 
presented for a bit-stuffed data source in Chapter V could be combined with the joint field- 
based maximum a posteriori detector (JFBMAPD), yielding even better FER performance. 

For signals that are received with low SNR, summing the soft outputs of a MAP 
detector may be used to increase the likelihood of successful demodulation. This is perhaps 
more relevant to single-channel low-entropy signals. 

In the area of parameter estimation, additional areas for further investigation in- 
clude joint carrier-phase and frequency tracking methods. This work assumes a source 
of channel state information (phase and frequency estimates). Determining this is non- 
trivial. In theory, these parameters could be treated as additional unknown parameters in 
the trellis. 

In the area of computation reduction, there are a few areas for further investigation. 
The optimum detector could be simplified by only searching back one symbol. This is 
similar to compressing the priors to one symbol. Looking back only one symbol is not 
optimal, but it would be interesting to compare the performance to Joint Field-based 
BCJR. There are also many variations on the BCJR algorithm that involve simplifying the 
computation [109]. 


149 


THIS PAGE INTENTIONALLY LEFT BLANK 


150 


APPENDIX A. LAURENT REPRESENTATION 


The ability to represent a signal as a linear combination of amplitude modulated 
pulses greatly simplifies receiver analysis and design. It allows one to build on the theory 
available for pulse amplitude modulated signals. 

Continuous-phase modulation (CPM) signals are in general nonlinear. CPM signals, 
in general, can not be represented as superposition of amplitude modulated pulses. Laurent 
[21] shows that for a finite length pulse response g(t), any CPM signal can be written as 


= gl-1 


a superposition of M amplitude modulated pulses, where L is the duration of the 


finite length pulse, g(t), divided by the symbol period. This is an extension of the well 
known representation of MSK as offset quadrature phase shift keying (OQPSK). We begin 


with the derivation of the Laurent representation and then follow with two examples. 


A. THE LAURENT REPRESENTATION 


The discussion below follows from examples in [21,23]. We start with the complex 


) (A.1) 


where / is the modulation index, 6 is an arbitrary constant phase term (assumed 0 here- 


envelope representation of a continuous phase modulated signal 


Oy + 2th ~ ang(t — nT) 


n 


z(t) = Aexp [ 








after), and A is a constant amplitude term (assumed 1 hereafter). The symbol a,, € {—1, 1} 
contains the information for the n“” bit. The function q(t) is the integrated finite duration 


pulse response 


0 ant 
q(t)=4 fo g(r)dr O<t<LT, (A.2) 
1/2 tL, 


for integer L and symbol duration T,. For GMSK, g(t) is a Gaussian pulse shape. For the 
case of AIS operating in the “high” setting with a time-bandwidth product of BT, = 0.4, 


L can safely be set to 4 without significant loss of energy. 


151 


Now we rewrite Equation (A.1) by extracting the phase contribution from all but 
the last L symbols. We then write the contribution of the last LZ symbols as the product 


of exponentials. 


z(t) = exp ta Ss o| x [lex (j27h |an_ig(t — (N — 1)T;)]) (A.3) 


i=0 


n=—Co 


The next step involves rewriting the complex exponential associated with the n” symbol 


as the sum of two terms with only the second term depending on the n“” symbol a, 


. 2 _sin(th — Ithq(t — nT;)) 
exp (j2rh lang(t Ts)]) sin(h) 








A.4 
Gin pomieanale — nT;)) ae 
Ren sin(7h) 
This can be done by recognizing we want the form 
sin(ah — B(t — nT)) + ysin(G(t — nT)) (A.5) 


and using trigonometric identities to solve for the coefficients y and 3. We then define a 


symmetric generalized phase pulse: 





2thq(t) a ew bd Ee 
= (A.6) 
2nh(1/2—q(t-—LT;)) LT, <¢t 
For notational convenience we define the function 
sin[W(t + nT,)| 
nm = = T; A. 
We can now substitute Equation (A.6) and (A.7) in Equation (A.4) to obtain 
exp (j [a, Y(t — nT;)]) = Szr_n(t) + exp(tha,)S_n(t) (A.8) 


We can now use Equation (A.8) in Equation (A.3) to get 


z(t) = exp [is Ss °| [[lSe+2-w(t) + exp(jahan_i)Si_n(t)] (A.9) 


i=0 


n=—CoO 
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The L product terms on the right side of Equation (A.9) can be expanded to a sum of 2” 
product terms, (actually 2’~! from symmetry). Following the notation in [22] we introduce 
the index K € {0,...,24-1 — 1}. We use the notation K(i) € {0,1} to represent the i” 


bit in the binary representation of K 
K=)_ K(i)2' (A.10) 


For example if L = 3 and K = 3, then K(0) =1, K(1) =1, and K(2) =0. 
We now group these terms defining what is commonly referred to as the Laurent 


pulses. 


Cx(t) = So(t) x I Si+K (iL (t) (A.11) 


Observe the Laurent pulses C'x(t) are formed from a product of shifted So(t) pulses. We 


then use the A(i) notation to define the complex phase coefficients 


N 
Aon = S- An (A.12) 
and 
Boa 
Ax,n = Aon — S- K(t)an— (A.13) 


i=1 
and express the original CPM signal as a sum of amplitude modulated pulses 
co «2h o1_] 


Zt) = S° S- exp(jthAx,w)Ck(t — nT) (A.14) 


N=-co Kk=0 
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B. EXAMPLES 


In this section we provide an example of the Laurent representation for MSK and 


an example of the Laurent representation for GMSK as used in AIS as defined in [1] using 


a reduced number of terms. 


1. Laurent Representation for MSK 


The Laurent representation for MSK leads to the well known interpretation of MSK 
as OQPSK with half sine waves as pulses. For MSK, h = 1/2, L = 1, and the pulse g(t) is 


given by 
a 062s 


g(t) = 4" 
0 otherwise 


The integrated pulse term q(t) is then 


0 t<0 


=o> USta i, 


1/2 otherwise 


The symmetric generalized phase pulse W(t) is given by 


a Apes ae ean 
WO) = 9 1 <t <2, 
0 otherwise 


We can use Equation A.7 to find So as 


sin || 0<t<27, 


0 otherwise 


Each of the above signals is shown in Figure A.1 
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(A.15) 


(A.16) 


(A.17) 


(A.18) 


























T. ui 
(a) g(t) (b) a(t) 
w(t) So(t) 
se QT. : Te QT. , 
(c) W(t) (d) So(t) 


Figure A.1: Pulse shapes for MSK 


For L = 1 the only Laurent pulse is the Cp term Co(t) = So(t) and the only complex 
phase coefficient is Ap,y. So for the case of MSK 


(oe) 


PA S- exp(jthAo,n)Co(t — nT) (A.19) 


N=—0oo 


which is identical to the OQPSK representation as a superposition of half sine wave pulses. 


2. Laurent Representation for AIS 

The Laurent representation for GMSK as used in AIS is significantly more involved 
than the representation for MSK. We will show that for GMSK with BT = 0.4, h = 0.5, 
and L = 4, although 16 terms are required for exact representation, using only the first 


one or two terms results in minimal error. 
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For GMSK, g(t) is a Gaussian pulse. 


1 = 
g(t) = rect(t/T,) * eye exp (sr (A.20) 


where * denotes convolution, ¢ = \/In(2)/27BT, [20], and rect(t) is defined as 


Jl it}<1/2 
rect(t) = (A.21) 
0 |t])>1/2 


The pulse response g(t) can be written as a difference of Q functions?, 


(sean $e) ay (sen tt) ren 





g(t) = oT. 


T;/In(2) Ts./In(2) 


where BT, is the time-bandwidth product of the pulse response g(t) and B is the 3dB 
bandwidth of the Gaussian filter. 
Typically g(t) will be truncated to length LZ depending on the value of BT, such 





that g(t) is zero outside the interval [0, LT,] and symmetric around LT,/2. 
The integrated pulse term q(t) is then 





0 tO 
at)= 4 fog(r)dr O<t<LT, (A.24) 
1/2 tL, 
The function W(t) is given by 
2rhq(t) O25 = LT; 
W(t) = ( —2xhq(t —2LT,) LT, <t < 2LT, (A.25) 
0 otherwise 
'The Q-function is defined as 
ae | 2 
— 7" /2 
Q(t) =| aa dr (A.22) 
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LT, : LT, 
(a) g(t) (b) a(t) 
vu) So(t) 
LT; 2LT, : LT; 2LT, : 
(c) W(t) (d) So(t) 


Figure A.2: Pulse shapes for GMSK 
We can use Equation A.7 to find So 


sin(W(t)) O<t< 2LT, 
0 otherwise 
Each of the above signals is shown in Figure A.2 
For L = 4 there will be M = 24"! = 8 Laurent pulses {C,(t),...,C7(t)} which can 
be found by evaluating Equation A.11. For illustration, the product terms to form the 8 
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Figure A.3: The first four pulses Co,C,,Co, and C3 for the Laurent representation with 
gaa ha )s, big 4 


pulses are shown below. 


CH= HKAO XO SSO 
Cy (t) =So(t) x Ss(t) x So(t) x S3(t) 
ORG) ASN OS) EMD ID 
C3(t) =So(t) x S5(t) x Se(t) x S3(t) (A.27) 
C4(t) =So(t) x S(t) x So(t) x Sr(t) 
C(t) =So(t) x S(t) x So(t) x S4(t) 
Co(t) =So(t) x S(t) x So(t) x Sr(t) 
C7(t) =So(t) x S5(t) x Se(t) x S7(t) 


Figure A.3 illustrates the first four Laurent pulses. Notice that the first four Lau- 
rent pulses {Co, Ci, C2, C3} contain 99.905%, 0.095%, 0.000014%, 0.000014% of the energy 


respectively. 
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The complex phase coefficients are calculated with Equation A.12 and A.13. For 


illustration, the 8 complex phase coefficients are shown below. 


N 
Ao,n = y An 


n=—oo 
N 
Ain _ y An — AN-1 
n=—oo 


N 
Aon = y An — AN-2 


n=—Co 


N 
A3,Nn = y An — AnN-1 — AN-2 


N 
(A.28) 
Agn = S- An — AN-3 


N 
As,N = y Qn — AN-1 — AN-3 


n=—OCoO 


N 
Aon = y An — An-2 — AN-3 


n=—Co 


N 
Az,N = y Qn — AN-1 — AN-2 — AnN-3 


n=—Co 


The complex part of z(t) formed using three methods: an integrated phase term, 
only the first Laurent pulse, and finally the first and second Laurent pulse. This is shown 
in Figure A.4. Notice that using only the first two pulses allows us to match that of the 
integrated phase term with a mean squared error (MSE) of only 2.3 x 107° %. Using only 
the first Laurent pulse results in a MSE of only 9.4 x 107? %. 
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Imaginary part 
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Signal using integrate phase term ——G—— 


= Signal using lst Laurent Pulse —--—-o--- 


Signal using lst and 2nd Laurent Pulse ----+---- 
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Figure A.4: GMSK signal generated with integrated phase term, first Laurent pulse and 
then the first and second Laurent pulse 
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APPENDIX B. BCJR ALGORITHM 


The MAP algorithm is an efficient algorithm for estimating the a posteriori proba- 
bilities of the state transitions of a Markov source. The MAP algorithm, also called BCJR 
after the authors Bahl, Cocke, Jelinek and Raviv was originally proposed in [83]. Consider 
a sequence of values a € A where A is a finite alphabet. The values a are then mapped 
one-to-one to symbols a and passed through a channel with finite impulse response. We 


can write the output of a whitened matched filter for this channel as 


pet 
ylr] = S0 afnjafe — n] + 2[n] (B.1) 
k=0 
where L is the length of the equivalent discrete-time channel h, and z is sampled bandlim- 
ited AWGN with variance o?. The MAP algorithm can be used to calculate the a posteriori 
probability Pr(a = a,|y), where a, is a specific element in A. 

Note: Throughout this dissertation, a has been used as above to refer to the trans- 
mitted symbol. For the remainder of this appendix a is used to represent an auxiliary 
probability function following the standard notation for BCJR. 

The MAP algorithm is similar to the Viterbi algorithm in that it takes into account 
an entire received sequence in making decisions. The key difference is the MAP algorithm 
also incorporates a priori information of the transmitted symbols and the outputs of the 
MAP algorithm are a posteriori probabilities for each symbol. Where Viterbi selects the 
maximum likelihood sequence given a received sequence, MAP generates the a posteriori 
probabilities for each symbol given a received sequence. 

The assumption for BCJR is the source can be written as a Markov process. Let 
A|n] be the state of the underlying Markov process at time n. The next state Aln + 1] is 
determined by the value of the input and the value of the previous state A[n]. If we know 
the starting state and all the inputs, we know all the state transitions. Similarly, if we 
know the state transitions, we know the inputs. 

Given a source that generates independent symbols with a set of known a priori 
probabilities passed through a memoryless channel, the MAP algorithm provides no ad- 
ditional benefit from a symbol-by-symbol MAP detector. The MAP algorithm allows one 
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to take advantage of the dependency among adjacent symbols. In the development of the 
BCJR algorithm below we draw on examples in [83], [109], and [108]. 

In our model, the “encoder” is a FIR filter. We encode the states such that the 
Markov property is maintained; i.e., knowledge of the state at time n+ 1 makes knowledge 
of the state at time n or previous times irrelevant. 

We start by rewriting the received vector y as the union of prior, current, and future 
observations. This allows us to exploit the Markov property to simplify the calculation of 


the posterior probabilities. 


Y = Yen UYn UYsn (B.2) 


where ye,, Yn; and y,,, indicate the prior, current and future observations respectively. 


We now define some auxiliary probability functions for notational convenience. 


On(p) = Pr(A[n] = p, Yen) (B.3) 


Equation (B.3) is the probability of being in state p at time n and all of the received values 


prior to n. 


Br+i(@) = Pr(ys,lAln + 1] = @) (B.4) 


Equation (B.4) is the probability of the received sequence of values received after n given 


the state at time n+ 1 is q. 


Yn(p, @) = Pr(Aln + 1] = @, yn |Aln] = p) (B.5) 


Equation (B.5) is the probability of the observation at time n and the state being q at time 
n+ 1 given a previous state at time n of p. 

Using Equation (B.2) and the identity Pr(A|B) = Pr(A, B)/ Pr(B) we can write 
the posterior probability as 


Pr(A[n] = p, A[n + 1] = gly) = PrQ[n] = p, Aln + 1] =, Yen Yn Ysn)/Prly) — (B.6) 
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We now factor (B.6) using the identity Pr(A, B) = Pr(A) Pr(B|A) to get 


Pr(A[n] = p, A[n + 1] = gly) =PrQ[n] = p, Aln + 1] = @ Yen: Yn) x 


(B.7) 
Pr(ysnlAln] = p, Alin + 1] = 9, Yens Yn)/ Pry) 
The first term in (B.7) can be factored using the same identity to get 
Pr(\[n] = p, A[n + 1] = gly) =PrO ln] = vp, yen) 
Pr(A[n + 1] = g, ynlA[n] = p, Yen) x (B.8) 


Pr(ys,|Aln] = p, Aln + 1] = @, Yen» Yn)/ Pr(y) 


The Markov property says that knowledge of the state at time n renders knowledge 
of previous states and previous inputs irrelevant. We can use this property on the second 
and third factors on the right hand side of (B.8) to get 


Pr(A[n] = p, A[n + 1] = gly) = PrQ [nr] = p, Yen) x 
Pr(A[r + 1] = 9, ylA[n] = p)x (B.9) 
Pr(ys,|A[n + 1] = q)/ Pr(y) 


We can now rewrite (B.9) in terms of the previously defined a, 3, and y. 


Pr(A[n] = p, A[n + 1] = gly) =an(p) x 
Yn(P; 4) x (B.10) 
Bn+i(q)/ Pr(y) 


There is a specific input associated with each state transition, but there may be 
multiple state transitions associated with each possible input. In other words, there is a 
many-to-one mapping of state transitions to inputs and a one-to-many mapping of inputs 
to state transitions. In order to find Pr(a = a,|y) we must sum over all the state transitions 


associated with the specific input a,. 


Pr(a =azly) = > an(p)yn(p, @)Gn41(9)/ Pry) (B.11) 


(p,q)€Sx 
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where S, is the set of transitions associated with the specific input a, 
We can remove the Pr(y) term by recognizing summing over all Pr(a = a,|y) for 


all a, € A will equal 1. 


y. Prieeacy) = (B.12) 


azEA 


We can now rewrite Equation (B.11) as 


Pr(a = ay|y) = ye (pg)eSe &n(P) InP, 2) Bn41(G)/ Pr(y) 
om ase A 2u(p.q)eSe On(P) In (Ps YGn+1(9)/ Pr(y) 





(B.13) 


= 2 (p,a)eSx On(P)Yn(P; 7) Bn+1 (4) 
ee DS pice @ On(p )In(P, DPn+1(9 ) 





where the denominator normalizes the probabilities to one. 
Both a@ and ( can be calculated recursively. The derivation below of the recursive 
calculation of a and ( follows [109]. We start by writing out the definition of a, (B.3) at 


the next time step. 


Ansi(q) = Pr(A[n + 1] = 9, Yens Yn) (B.14) 
the marginal distribution is computed from the joint, Pr(A) = }) pe, Pr(A, B) 


Q-1 
On+1(q) = Pr(A[n ao 1] = 4, Yen Uns A[n] = p) (B.15) 


p=0 
where @ is the number of possible states at time n. This is a sum over all possible previous 
states. 
We now factor the right hand side of Equation B.15 using the identity Pr(A, B) = 
Pr(A) Pr(B|A) to get 


Onti(q) = )»_ Pr(A[n] = p, Yen) Pr(Aln + 1] = 4, mnlAln] = PD, Yen) (B.16) 
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From the Markov property if we know the state at time n, knowledge of values of 


y before time n are irrelevant. 


Q-1 
An+i(q) = D>, Pr(Aln] = Pp, Yen) Pr(Aln + Y= @, yn|Aln] = p) 
-_e-eoeo-O eo Se 
p=0 
an(p) "In (P.9) (B.17) 
Q-1 
p=0 


An almost identical method can be used to find a backward recursion for determining 


G resulting in 


Q- 
= TP An + 1] = 4¢, yn|Aln] = p) Pr(ys,, |Aln + 1] = @) 
AU gh 5 Pe i 
= Yn (p,q) Bn+1(9) (B.18) 
Q-1 
In (p q) Bn+1( q) 
q=0 


where @ is the number of states in the next state at time n+ 1. 
The term y represents the probability of transition from p to q at time n with the 


observation y,. We can factor 


Yn(p, gd) =Pr(A[n + 1] = 4, ynlAln] = p) 
=Pr(ynlAln + 1] = @,Aln] =p) PrAln+ 1 =alaln]=p) (8-19) 


a priori probabilities 


The second factor is the a priori probability Pr(a, = a). The first factor can be rewritten 
as Pr(Yn|Spq) Where sp is the output of the channel uniquely associated with a state 


transition from p to gq. For an AWGN channel we can write this as 


202 


1 —llyn = Spall? 
Privy |siq):= noi oP Ga (B.20) 


where o? is the variance of the noise term from Equation B.1. 
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Table B.1: MAP Decoding Algorithm 
Qo <= initial a 
Bn <= initial 
for n = 0 to N —2 do 
for q=0 to Q-—1do 


Onti(g) = Ng Dyno On(p)yn(P, 9) {Forward} 
end for 


end for 
for n = N —2 to 0 do 
for p= 0 to Q—1 do 


Br(p) = Np ary Yn(P, @)Bn+1 (9) {Backward} 
end for 


end for 
for n =0 to N —2 do 
for all a, € Ado 


Pr(a = azly) = Naz Xi(pqes, On(P) In(P; )Fn+1(4)/ Pr(y) {a posteriori} 
end for 


end for 


For numerical stability reasons, it is common to normalize both @ and (@ at ev- 
ery stage. The effect of the normalization cancels out when calculating the a posteriori 
probabilities with Equation (B.13) 

A summary of the MAP algorithm for calculating the a posteriori probabilities is 
given in Table B.1. The N, factors in front of the summations are normalizations for each 
step such that the as and (@s sum to one for each stage. 

In summary the MAP algorithm is an efficient method for calculating the a pos- 
teriori probabilities. Observe from Algorithm B.1 that the calculation of the a posteriori 


probabilities involves three passes, each pass with roughly N x Q? operations. 
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